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1.  Progress  of  Theoretical  Work 

In  the  AFOSR  grant,  the  complex  interaction  between  noise,  stability 
and  nonlinear  dynamics  inherent  in  mechanical  systems  was  examined 
and  a  consistent  method  for  analyzing  stochastic  nonlinear  systems  was 
developed.  This  has  led  to  a  new  qualitative  understanding  of  the  physical 
phenomena  encountered  in  nonlinear  mechanical  and  structural  systems 
under  harmonic  and  stochastic  excitations.  The  mathematical  techniques 
developed  under  the  current  AFOSR  grant  will  allow  engineers  and 
scientists  to  predict  possible  motion  instabilities  in  such  systems  as  rotor 
blade  dynamics  in  forward  flight  in  a  region  of  high  atmospheric 
turbulence  as  well  as  rotating  shafts  subjected  to  harmonic,  stochastic  and 
combined  harmonic  and  stochastic  excitations. 

1.1  Nonlinear  Deterministic  Systems 

Recently,  considerable  effort  has  been  directed  toward  obtaining  a 
better  understanding  of  nonlinear  behavior  and  instability  mechanisms  of 
rotating  shafts,  a  fundamental  component  of  many  mechanical  systems. 
Toward  this  end,  an  analytical  method  based  on  both  Hamiltonian  and  non- 
Hamiltonian  frameworks  is  being  developed  by  the  PI  and  one  of  his 
graduate  students.  The  case  in  which  a  harmonic  axial  excitation  is 
applied  has  been  studied  in  detail  by  the  PI  and  his  graduate  student  [1],  It 
was  shown  that  in  addition  to  the  simple  and  Hopf  bifurcations  in  the 
presence  of  subharmonic  and  combination  resonances,  respectively,  such 
systems  can  exhibit  a  generalized  Hopf  bifurcation  with  1:1  resonance 
when  the  damping  is  very  small.  In  connection  to  this,  the  effects  of 
periodic  parametric  excitations  on  systems  exhibiting  Hopf  bifurcations 
with  1:1  resonance  are  also  being  investigated  [2,  3].  Both  semisimple  and 
non-semisimple  cases  are  of  interest,  and  the  normal  forms  are  calculated 
when  the  forcing  frequency  is  near  twice  the  natural  frequency. 
Furthermore,  since  symmetric  properties  are  fundamental  features  of 
many  engineering  systems,  the  PI  is  currently  studying  global  bifurcations 
with  rotational  and  reflective  symmetries,  as  well  as  the  effects  of 
imperfections  that  destroy  this  symmetry. 
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The  techniques  developed  under  the  current  AFOSR  grant  can  be 
applied  to  a  wide  variety  of  realistic  problems  in  nonlinear  structural 
dynamics.  As  mentioned  above,  the  results  have  already  been  applied  to 
study  the  stability  and  bifurcation  behavior  of  gyroscopic  systems  (i.e.  the 
rotating  shaft)  under  combined  harmonic  and  stochastic  excitations  and 
the  dynamics  of  rotor  blades  in  forward  flight.  Further  applications 
include  the  analysis  of  propellant  lines  conveying  pulsating  (possibly 
turbulent)  fluid  flow.  The  PI  has  introduced  novel  mathematical  and 
computational  methods  to  study  the  nonlinear  behavior  of  both  gyroscopic 
and  nonconservative  deterministic  systems  [4,  5].  Over  the  next  few  years, 
the  PI  will  develop  both  local  and  global  techniques  to  investigate  various  co¬ 
dimension  two  and  three  bifurcations  under  time  periodic  perturbations. 
Special  emphasis  will  be  given  to  the  study  of  global  chaotic  phenomena  in 
these  systems. 

1.2.  Nonlinear  Stochastic  Systems 

In  many  situations,  parametric  or  external  excitations  cannot 
always  be  adequately  modelled  by  deterministic  time  functions  alone.  They 
fluctuate  randomly  over  a  wide  band  of  frequencies  and  have  to  be 
considered  as  stochastic  functions  of  time  defined  only  in  probabilistic 
terms.  Since  the  effects  of  stochastic  perturbations  are  of  greatest 
importance  near  bifurcation  points  in  any  dynamical  system,  a  portion  of 
this  research  focuses  on  noise  induced  transitions  near  such  points.  The 
results  obtained  have  great  impact  on  such  engineering  problems  as 
aircraft  at  high  angles  of  attack  [61,  panels  under  gas  flow  with  both 
turbulent  boundary  layers  and  random  axial  loads  [7,8],  rotating  systems 
under  pulsating  axial  loads  [9]  and  propellant  lines  conveying  turbulent 
fluid  flow  [10]. 

The  PI  has  developed  a  method  called  "stochastic  normal  forms” 
which  replaces  the  original  nonlinear  system  by  an  "equivalent”  (in  the 
stochastic  sense)  system  of  lower  dimension  [11,12],  A  consistent  method 
for  analyzing  stochastic  nonlinear  systems  has  been  developed  through  this 
investigation.  These  new  techniques  [13,  14,  15,  16],  based  on  the  concept  of 
Lyapunov  exponents  and  multiplicative  ergodic  theory,  provide  insight  into 
the  effects  of  harmonic  and  stochastic  excitations  on  the  response  of 
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nonlinear  dynamical  systems.  Stability  boundaries  and  bifurcation 
scenarios  are  then  determined  for  such  systems  as  helicopter  rotor  blades 
in  forward  flight  [17,  18j  and  rotating  shafts  subjected  to  combined 
harmonic  and  stochastic  excitations  [9].  Results  obtained  in  this  study 
explain  how  small  amplitude  periodic  or  stochastic  fluctuations  in  the 
parameters  of  a  system  or  its  environment  can  have  a  marked  effect  on  the 
dynamics  of  physically  realistic  nonlinear  systems.  The  current  research 
has  led  to  fundamental  contributions  in  understanding  co-dimension  two 
stochastic  bifurcations  such  as  the  interaction  of  Hopf  and  pitchfork 
bifurcations  under  stochastic  excitations  [19,  20]. 

In  the  prior  AFOSR  support,  the  results  were  obtained  for  white  or 
nearly  white  noise  cases.  These  assumptions,  however,  are  idealistic.  The 
PI  is  currently  examining  the  the  cases  of  realistic  colored  noise  and 
simultaneous  harmonic  and  stochastic  excitations  which  are  often 
encountered  in  nonlinear  mechanical  and  structural  systems. 

2.  Progress  of  Experimental  Work 

The  objective  of  the  experimental  research  is  the  development  of 
laboratory  facilities  to  conduct  a  series  of  tests  designed  to  verify  the 
analytical  results  obtained  by  the  PI  under  his  current  AFOSR  grant  and 
other  previous  support. 

2. 1  Mechanical  Problem 

One  of  the  most  fundamental  components  of  a  mechanical  system  is 
a  rotating  shaft.  It  is,  therefore,  not  surprising  that  through  the  years 
considerable  effort  has  been  directed  toward  obtaining  a  better 
understanding  of  such  mechanisms.  The  dynamics  and  response  of 
rotating  and  gyroscopic  systems  have  been  studied  extensively  in  the 
literature.  The  equations  describing  the  transverse  motion  of  a  continuous 
uniform  elastic  shaft  of  asymmetrical  cross-section  mounted  in  a  rigid 
bearing  and  rotating  with  angular  velocity  £2(t)  about  the  horizontal  center 
line  of  the  bearings  are 
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EIU +  P(t)  ~  +  m ~  +  D^-  -  2mD(t)|^  -  mQ(t)2u  =  fi{u,v) 

dz4  9z2  8t2  *  ox 

i-.t  .  d\  9 ~ 9v  _  ,9u  ,  .2  - ,  , 

Ely  —  +  P(t) ~  +  m  — +  D-r-  -  2mfi(tW~ -  mQ(t)  v  =  ftu.v) 

dz4  9z2  8t2  ox  8t 

where  P(t)  may  be  harmonic,  stochastic  or  a  combination  of  both.  The  shaft 
is  of  length  L,  mass  per  unit  length  m,  and  flexural  rigidities  EIU  and  EIV 
parallel  to  directions  Ou  and  Ov,  respectively.  The  terms  fj(u,v),  i=l,2, 
consist  of  all  nonlinearities  due  to  nonlinear  strain-displacement  relations 
and  the  nonlinear  damping  model.  In  addition,  fj(u,v)  may  also  contain  the 
effects  of  mass  imbalance  and  bearing  forces,  including  impacts. 

2.2  Proposed  Experimental  Investigation  of  the  Rotating  Shaft 

One  of  the  main  objectives  of  the  proposed  research  is  to  conduct  a 
series  of  experiments  on  a  rotating  shaft.  An  important  application  of  this 
work  is  the  study  of  aircraft  gas  turbine  engines.  Both  rotor  imbalance  and 
impact  due  to  bearing  clearance  are  inherent  in  such  systems.  Thus,  the 
experimental  set-up  will  consist  of  an  elastic  shaft  supported  by  bearings 
with  clearance  to  depict  realistic  operating  conditions  in  a  turbine  engine. 
This  clearance  is  modelled  by  a  bilinear  spring.  It  is  intended  to  measure 
the  time  response  of  a  rotating  shaft  under  two  separate  operating 
conditions.  Specifically,  the  transverse  vibrations  will  be  monitored  as  the 
shaft  is  subjected  to  (1)  time  varying  axial  thrust  at  various  constant 
rotation  rates  and  (2)  various  static  axial  loads  with  time-dependent 
angular  velocity. 

The  experimental  problem  can  be  divided  into  three  categories: 
design  and  construction  of  the  shaft  and  test  rig,  integration  of  the  shaker 
(external  excitation)  and  motor  (rotation),  and  implementation  of  sensors 
and  time  series  analysis.  Detailed  descriptions  of  each  of  these  components 
are  given  in  the  following  four  subsections. 

2.2.1  Design  of  Shaft  and  Test  Rig 

The  shaft  to  be  studied  will  be  constructed  from  hard  naval  brass. 
The  published  material  properties  of  hardened  yellow  brass  are  shown  in 
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Table  1  along  with  the  shaft  characteristics  of  interest.  The  shaft  has  a 
clamped  length  of  L  =  25.1  in.  and  a  diameter  of  d  =  0.25  in. 

As  shown  in  Table  1,  the  critical  static  buckling  load  of  the  shaft  is 
Pcr  =  190  lbf  and  the  natural  angular  frequency  is  fn  =  50  Hz.  Assuming  a 
safety  factor  of  1.5  (amax  =  Cyjeid/1.5),  the  maximum  allowable  strain  will  be 
emax  =  2500  pin/in  which  corresponds  to  a  centerline  deflection  of  a  =  0.3  in. 


Table  I:  Physical  Characteristics  of  Shaft 

Shaft  Property 

Characteristic  Value 

Material 

hard  naval  brass 

E 

15.9  x  106  psi 

M- 

7.92  x  HH  lbf-s2/in4 

sy 

60  ksi 

L 

25.1  in 

d 

0.25  in 

fn 

50  Hz 

Per 

191  lbf 

amax 

0.3  in 

em«x 

2500  pin/in 

(Ax)max 

0.018  in 

The  bearing  assemblies  are  designed  to  closely  approximate 
clamped-clamped  boundary  conditions.  With  the  use  of  super-precision 
bearings,  the  maximum  angular  deflection  of  the  shaft  at  the  bearings  is 
calculated  to  be  0.013  degrees.  Thus,  for  the  first  0.013  degrees  of  deflection 
at  the  bearings,  the  model  behaves  as  if  it  is  a  pinned-pinned  shaft. 
Assuming  pinned-pinned  boundary  conditions,  the  maximum  center  line 
deflection  corresponding  to  this  angular  deflection  at  the  ends  can  be 
calculated  to  be  0.002  inches.  This  deflection  corresponds  to  a  strain  gage 
reading  of  approximately  1.9  pin/in  which  is  negligible.  Thus,  we  can 
assume  clamped-clamped  boundary  conditions  throughout. 

The  brass  shaft  is  held  in  the  motor  end  bearing  hub  using  a  slit 
clamp  and  bolts  and  four  set  screws.  The  large  steel  bearing  blocks  and 
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large  thrust  bearing  take  the  full  load  transmitted  by  the  shaft.  The  bearing 
hub  has  eight  port  holes  tapped  through  it  to  allow  for  passage  of  the  strain 
gage  leads  from  the  shaft  to  the  slip  ring,  which  is  mounted  on  the  bearing 
hub.  The  bearing  hub  is  then  passed  through  a  radial  bearing  and  coupled 
to  the  flywheel  and  drive  motor.  Also  mounted  on  the  bearing  hub  is  a 
timing  belt  pulley  used  to  drive  the  auxiliary  shaft.  The  auxiliary  shaft  and 
timing  belts  are  used  to  synchronize  the  rotation  rate  of  the  shaker  end 
bearing  hub  with  that  of  the  shaft  to  minimize  torque  in  the  experimental 
model.  Therefore,  the  shaft  will  experience  only  axial  movement  with 
respect  to  the  linear/rotary  bearings,  insuring  that  they  do  not  bind. 

The  shaker-end  bearing  hub  is  mounted  in  two  bearings,  one  of 
which  is  a  thrust  bearing  used  to  provide  increased  axial  rigidity.  Three 
0.375  inch  linear/rotary  bearings  are  mounted  inside  the  bearing  hub  to 
allow  axial  movement  of  the  shaft.  The  hardened  stainless  steel  sleeve  on 
which  the  linear/rotary  bearings  act  is  attached  to  the  shaft  using  a  shrink 
fit.  The  sleeve  acts  as  the  inner  race  for  the  ball  bearings  and  also 
increases  the  stiffness  of  the  shaft  through  the  bearing. 

A  small  thrust  bearing  is  used  to  transmit  the  axial  loading  to  the 
rotating  shaft.  A  load  cell  is  mounted  on  the  thrust  bearing  to  measure  the 
static  and  dynamic  axial  loads  which  are  to  be  applied.  Mounted  directly 
above  the  load  cell  is  a  cross  bar  which  transmits  the  static  and  dynamic 
axial  loading  to  the  shaft.  The  static  load  is  applied  using  two  small 
pneumatic  air  cylinders  which  are  pressurized  with  a  one  gallon  air  tank. 
An  air  tank  pressure  of  100  psia  will  apply  a  static  load  of  260  ibf.  The 
dynamic  axial  load  is  applied  to  the  shaft  using  a  100  lbf  shaker,  which  acts 
through  the  center  of  the  cross  bar. 

2.2.2  Integration  of  Motor  and  Shaker 

The  motor  currently  under  consideration  for  these  experiments  is  the 
MTS  1000-Watt  pulse-width  modulated  servo  motor  and  controller.  This 
motor/controller  system  provides  high  performance,  high  accuracy 
operation  necessary  for  the  proposed  experimental  project.  In  the  case  of  a 
constant  rotation  rate,  i.e.  fl(t)=ft0,  such  a  high  torque  motor  is  not 
necessary.  However,  in  the  second  phase  of  the  proposed  experimental 
program,  the  parametric  excitation  is  given  in  the  form  of  a  time  dependent 


8 


rotation  ra  i.e.  Q(t)  =  ft0  +  Gi(t).  The  high  torque  required  to  insure 
precise  control  of  the  time-varying  component  of  this  excitation  necessitates 
a  high  performance  motor. 

The  dynamic  axial  loading  applied  to  the  shaft  is  generated  by  a 
Bruel  &  Kjaer  445  N  (100  lbO  shaker  mounted  as  shown  in  Figures  3  and  5. 
The  shaker  is  used  to  apply  only  the  dynamic  load  since  use  of  this  system 
to  apply  static  loads  may  cause  the  shaker  to  overheat.  The  Bruel  &  Kjaer 
shaker  system  was  chosen  since  it  provides  the  most  precise  frequency 
command  following  which  is  essential  in  these  experiments.  The 
command  signal  is  generated  by  a  Tektronix  2630  FFT  Analyzer/Signal 
Generator.  This  signal  generator  is  capable  of  producing  dc,  periodic, 
random  or  any  combination  of  these  waveforms. 

The  amplitude  of  the  static  load  and  the  amplitude  and  frequency  of 
the  dynamic  load  will  be  measured  by  a  200  lbf  S-type  load  cell  mounted  on 
the  thrust  bearing  at  the  shaker  end  of  the  shaft.  The  signal  from  the  load 
cell  is  passed  through  a  Wagner  model  #460  Signal  Conditioner/Amplifier. 
Currently,  the  shaker  system  is  run  open-loop;  the  signal  from  the 
amplifier  is  input  to  the  Tektronix  2630  analyzer  where  it  can  be  compared 
to  the  reference  input.  In  this  configuration,  no  corrective  action  is  taken  if 
the  excitation  amplitude  does  not  correspond  exactly  to  the  desired 
waveform  (for  periodic  inputs)  or  to  the  desired  spectrum  (for  random 
inputs). 

2.2.3  Data  Acquisition  and  Time  Series  Analysis 

A  schematic  diagram  of  the  experimental  set-up  is  given  in  Figure  9. 
The  deflection  of  the  shaft  in  the  principle  directions,  i.e.  the  rotating  u-  and 
v-coordinates,  will  be  measured  by  strain  gages.  A  total  of  four  foil  gages 
are  mounted  on  the  shaft,  one  every  90°,  and  arranged  in  two  half-bridge 
configurations.  The  strain  gage  signals  are  passed  through  the  Fabricast 
model  #1401001  slip  ring  to  a  Measurements  Group  #2311  Signal 
Conditioning  Amplifier.  The  conditioned  signals  are  then  input  to  the 
Tektronix  analyzer.  The  analyzer,  run  from  a  Gateway  486  PC, 
implements  the  following  standard  signal  analysis  functions: 
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♦  time  domain  waveform  and  orbit  (x  -  y)  plots; 

♦  averaged  power  spectral  density  and  cross  spectral  density  functions; 

♦  transfer  (frequency  response)  function  and  FFT; 

♦  coherence  function; 

♦  waveform  averaging; 

♦  auto-  and  cross-correlation; 

♦  impulse  response  function; 

These  standard  functions  are  indeed  the  most  frequently  applied 
techniques  for  understanding  experimental  data.  Such  analysis  is  useful 
for  obtaining  the  frequency  components  and  power  distribution  as  a 
function  of  frequency.  However,  in  order  to  gain  a  detailed  description  of 
the  nonlinear  dynamics,  it  is  imperative  to  utilize  nonstandard  time  series 
analysis.  Lyapunov  exponents,  wavelet  transforms  and  fractal  dimensions 
are  examples  of  such  nonstandard  techniques. 

3.  Current  Research 

3.1  Specific  Objectives 

The  specific  objective  of  the  current  work  is  to  examine,  both 
theoretically  and  experimentally,  the  nonlinear  response  of  gyroscopic 
systems  under  stochastic  and  harmonic  excitation.  This  study  will  include 
the  following: 

♦  examination  of  the  combined  effects  of  mass  imbalance,  asymmetry 
and  realistic  boundary  conditions  on  the  nonlinear  response  of 
rotating  systems; 

♦  development  of  a  nonlinear  control  technique  designed  to  suppress 
unwanted  vibrations  and  chaotic  motions  in  gyroscopic  systems; 

♦  experimental  verification  of  the  theories  developed. 

The  current  research  will  investigate  the  effects  of  both  parametric 
and  internal  resonances  on  the  local  and  global  dynamics  of  realistic 
gyroscopic  systems.  In  most  realistic  cases,  the  mechanical  system  is 
subject  to  stochastic  as  well  as  deterministic  excitations.  The  stability 
boundaries  for  such  cases  must  also  be  obtained.  Toward  this  end,  the  PI 
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and  his  graduate  students  are  extending  the  techniques  developed  in  the 
study  of  stochastic  dynamical  systems  by  removing  some  of  the  restrictive 
conditions.  The  goal  of  this  effort  is  the  development  of  a  systematic 
approach  to  nonlinear  stochastic  problems.  The  results  obtained  will 
identify  the  mechanisms  which  give  rise  to  instability  and  unwanted 
vibrations  under  complex  loading  conditions.  In  addition  to  understanding 
the  mechanisms  of  instability,  it  is  beneficial  to  be  able  to  suppress 
unwanted  component  oscillations.  A  general  approximate  nonlinear 
control  methodology  designed  to  suppress  undesired  motions  in 
mechanical  system  such  as  the  rotating  shaft  is  also  under  development  as 
part  of  the  current  research.  In  the  above-mentioned  studies,  benchmark 
experiments  are  necessary  to  validate  the  current  theories  for  the  limited 
classes  of  systems  to  which  they  apply  and  to  guide  the  development  of  more 
general  theories  in  the  areas  of  nonlinear  dynamics  and  control. 

3^2  Scope  of  the  Current  Theoretical  Work 

The  research  in  the  previous  AFOSR  grant  has  produced  important 
analytical  results  for  lower  dimensional  systems.  However,  little  analysis 
and  few  results  auc  available  for  large  realistic  multidegree-of-freedom 
(n>4)  mechanical  systems  due  to  their  complexity.  Furthermore,  it  is 
imperative  that  experimental  results  be  provided  for  comparison  to  the 
theoretical  models.  This  project  addresses  both  the  analytical  and 
experimental  aspects  of  the  Pi's  current  and  planned  research. 

In  the  course  of  the  current  study,  several  questions,  which  have 
hitherto  not  been  studied  globally,  will  be  addressed: 

♦  What  are  the  effects  of  mass  imbalance,  asymmetry  and  motion 
constraints  on  the  nonlinear  response  of  gyroscopic  systems? 

♦  What  are  the  mechanisms  of  global  bifurcations  in  such  systems? 

♦  Why  is  noise  beneficial  in  some  nonlinear  systems  but  harmful  in 
others? 

♦  What  additional  effects  are  caused  by  the  presence  of  harmonic 
excitation? 

♦  What  is  the  effect  of  noise  on  the  stable  and  unstable  periodic 
motions? 


Answers  to  these  questions  will  deepen  our  understanding  of  the  dynamics 
of  mechanical  systems  in  real  life  situations  and  will  aid  in  the  design  of 
more  efficient  dynamical  systems. 

The  Pi's  current  research  is  also  concerned  with  the  development  of 
&  general  control  algorithm  for  nonlinear  multidegree-of-freedom 
mechanical  systems.  A  variety  of  power  generating  components  become 
unstable  due  to  parametric  excitations  which  may  lead  to  catastrophic 
failures.  Thus,  it  is  important  to  study  the  instability  mechanisms  in  detaii 
and  to  control  undesirable  component  motions.  This  work  will  result  in  a 
systematic  method  of  suppressing  unwanted  vibrations  and  chaotic 
motions  in  rotating  shafts  subjected  to  motion  constraints.  The  nonlinear 
control  designs  will  then  be  implemented  in  the  laboratory. 

3JJ  Scope  of  the  Current  Experimental  Work 

At  this  stage  of  our  investigation,  it  is  imperative  to  verify  the 
theoretical  results  of  the  ongoing  APOSR  grant  through  a  series  of 
experiments.  To  this  end,  the  PI  will  conduct  experiments  with  models 
whose  essential  dynamics  correspond  to  the  class  of  systems  for  which  the 
theory  was  developed.  These  experiments  are  intended  to  validate  the  new 
theories  and/or  to  identify  the  essence  of  physical  phenomena  that  must  be 
modelled. 

Unlike  the  area  of  deterministic  dynamics,  there  are  very  few 
experimental  studies  on  nonlinear  stochastic  dynamics  available  for 
comparison.  Therefore,  the  current  experimental  work  will  be  valuable  in 
providing  insight  into  the  stochastic  dynamics  of  actual  mechanical 
systems  and  in  bridging  the  gap  between  the  existing  theoretical  analysis 
and  physical  observation. 

The  numerical  and  analytical  results  will  be  compared  with  the 
folio  ing  estimated  results  obtained  from  the  experiments: 

♦  mean  squares  and  power  spectra  of  responses; 

♦  auto-  and  cross-correlation  of  the  response  coordinates; 

♦  probability  density  functions  of  the  responses; 

♦  Lyapunov  exponents  and  fractal  dimensions. 


The  experiments  will,  in  turn,  guide  the  developments  and 
refinements  of  the  theories  developed  to  incorporate  any  new  phenomena 
observed.  In  addition,  there  is  a  primary  need  to  increase  the  experimental 
skills  of  dynamics  and  control  graduate  students  and  the  current  program 
will  establish  a  modern  dynamics  and  control  laboratory  at  the  graduate 
level. 
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Normal  Form  for  Generalized  Hopf  Bifurcation 
with  Non-semisimple  1:1  Resonance1- 


N.  Sri  Namachchivaya  and  Monica  M.  Doyle 
Department  of  Aeronautical  and  Astronaut  ical  Engineering 
University  of  Illinois  at  Urbana-Champaign 
Urbana,  IL,  61810 

and 


William  F.  Langford  and  Nolan  W.  Evans 
Department  of  Mathematics  and  Statistics 
University  of  Guelph 
Guelph  Ontario,  Canada,  NIG  2W1 


The  primary  result  of  this  research  is  the  derivation  of  an  explicit  formula  for 
the  Poincare-Birkhoff  normal  form  of  the  generalized  Hopf  bifurcation  with  non- 
semisimple  1:1  resonance.  The  classical  nonuniqueness  of  the  normal  form  is 
resolved  by  the  choice  of  complementary  space  which  yields  a  unique  equivariant 
normal  form.  The  4  leading  complex  constants  in  the  normal  form  are  calculated 
in  terms  of  the  original  coefficients  of  both  the  quadratic  and  cubic  nonlinearities 
by  two  different  algorithms.  In  addition,  the  universal  unfolding  of  the  degenerate 
linear  operator  is  explicitly  determined.  The  dominant  normal  forms  are  then 
obtained  by  rescaling  the  variables.  Finally,  the  methods  of  averaging  and  normal 
forms  are  compared.  It  is  shown  that  the  dominant  terms  of  the  equivariant 
normal  form  are,  indeed,  the  same  as  those  of  the  averaged  equations  with  a 
particular  choice  for  the  constant  of  integration. 
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NSERC  of  Canada. 
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I.  Introduction 


When  a  multidegree-of-freedom  dynamical  system  undergoes  a  bifurcation,  it 
usually  does  so  in  only  a  few  degrees  of  freedom.  One  simple  example  is  the  buckling 
of  a  column.  If  p  and  p<;  represent  the  axial  and  Euler  loads  of  a  column,  respectively, 
then,  as  p  is  varied  in  the  vicinity  of  pc,  the  temporal  evolution  of  the  motion  is 
dominated  by  the  critical  mode  which,  in  the  first  approximation,  is  governed  by 
x  *  (p  -  ^  x  +  ax3.  A  more  complicated  situation  arises  when  several  control 
parameters  p  are  varied  in  such  a  way  that  several  modes  become  marginally 

unstable  simultaneously.  In  the  latter  case,  the  system  is  said  to  undergo  a  multiple 
bifurcation.  The  simplest  and  smallest  number  of  equations  which  capture  the 
essential  dynamics  of  the  original  system  in  the  vicinity  of  Pc  are  said  to  be  in  the 

normal  form.  The  theory  of  normal  forms  is  an  important  analytical  tool  for 
investigating  the  qualitative  behavior  of  nonlinear  dynamical  systems. 

The  idea  of  normal  forms  for  nonlinear  systems  dates  back  as  far  as  Euler; 
however,  Poincare  [16]  and  Birkhoff  [3]  were  the  first  to  bring  forth  the  theory  in  a 
more  definite  form.  Poincare  [16]  considered  the  problem  of  reducing  a  system  of 
nonlinear  differential  equations  to  a  system  of  linear  ones;  namely, 

dj-  *  Ax  +  f(x)  to  =  Ay,  xe  Rn,  ye  Rn.  Q) 

at  at 

The  formal  solution  of  this  problem  entails  finding  near-identity  coordinate 
transformations,  x  =  y  +  d>(y),  which  eliminate  the  analytic  expressions  of  the 
nonlinear  terms.  It  has  been  shown  that  such  a  formal  solution  exists  provided  the 
above  system  is  hyperbolic  and  the  eigenvalues  Xj  of  the  diagonalizable  matrix  A 
satisfy  the  nonresonance  condition 

Xi  *  X  mih  fori  =  1,2,...,  n  ,  ImlsX™1!^2  (2) 
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where  m  is  a  vector  of  integers  m  =  (mi ,  m2 ,...,  nan)  with  m/ >  0.  Furthermore,  it 
was  proven  that  if,  in  addition  to  the  above  results,  the  eigenvalues  lie  strictly  to  one 
side  of  a  line  separating  them  from  zero  in  the  complex  plane,  then  the  formal  series 
4>(y)  is  convergent. 

If  the  system  is  nonhyperbolic  or  condition  (2)  is  violated,  the  analytic 
expressions  of  the  nonlinear  terms  cannot  be  completely  eliminated  via  a  nonlinear 
change  of  coordinates.  The  remaining  terms  comprise  the  normal  form  of  the  system 
of  equations  given  by  (1).  The  normal  form  is  dictated  by  the  nature  of  the  linear 
operator  A.  Thus,  the  nonlinear  system  in  Eq.  (1)  can  be  reduced  to 

*  Ay  +  g  (y)  ,  ye  Rn  (3) 

where  g  is  simpler  than  f.  Such  reductions  have  been  widely  used  to  study 
deterministic  autonomous  and  nonautonomous  systems  (see  Arnold  [1]). 

In  bifurcation  problems,  the  eigenvalues  of  the  linear  operator  A  are  composed 
of  two  sets,  one  on  the  imaginary  axis  and  the  other  with  strictly  negative  real  parts. 
The  linear  vector  space  E  associated  with  A  can  also  be  divided  accordingly  as 
E  =  Ec  ®  Eg  such  that  xc  e  Ec  and  x8  €  Es  with  x  =  xc  +  xB.  There  are  two  approaches 
to  obtaining  normal  forms  for  deterministic  systems.  In  the  first,  as  shown  in 
Guckenheimer  and  Holmes  [12],  one  first  computes  the  lower  dimensional  center 
manifold  onto  which  the  dynamics  settle  for  large  times.  The  dynamical  system 
defined  on  the  center  manifold  is  then  transformed  to  the  normal  form  through  a 
nonlinear  change  of  coordinates.  In  the  second  method,  one  systematically  expands 
the  original  vector  field  in  powers  of  amplitudes  of  the  critical  modes  to  yield  both  the 
normal  form  and  center  manifold,  simultaneously,  as  shown  by  Elphick  et  al.  [7],  The 
approach  adopted  in  this  paper  for  the  computation  of  the  normal  form  assumes  that 
the  center  manifold  theorem  has  been  applied  to  the  original  system  and  is  based 
heavily  on  the  work  of  Elphick  et  al.  [7]. 
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The  aim  of  this  paper  is  two-fold:  first,  to  present  an  explicit  formula  for  the 
normal  form  of  a  generalized  Hopf  bifurcation  with  non-semisimple  1:1  resonance 
and,  second,  to  compare  the  results  with  those  obtained  via  the  method  of  averaging. 
The  results  for  the  corresponding  semisimple  case  were  obtained  by  Bajaj  and  Sethna 
[2]  using  center  manifold  theory  and  the  method  of  integral  averaging. 

Recently,  the  normal  form  for  a  generalized  Hopf  bifurcation  was  expressed  as 
a  4-dimensional  real  system  by  Cushman  and  Sanders  [5]  and  as  a  2-dimensional 
complex  system  by  Elphick  et  al.  [7]  and  Iooss  and  Adelmeyer  [18].  loose  et  al.  [19] 
employed  the  2-dimensional  normal  form  given  in  [7]  to  examine  the  steady 
bifurcating  solutions  in  nonlinear  hydrodynamic  stability  problems.  However,  there 
are  no  explicit  formulas  relating  the  coefficients  of  the  original  system  to  those  of  the 
normal  form.  This  paper  presents  explicit  formulas  for  the  4  leading  constants  in  the 
complex  normal  form  in  terms  of  coefficients  of  the  original  nonlinear  system  v.  :th 
both  quadratic  and  cubic  nonlinearities.  The  complex  normal  form  presented  by 
Elphick  et  al.  [7]  has  recently  been  analyzed  by  van  Gils  et  al.  [11],  It  was  shown 
that  this  co-dimension  3  bifurcation  problem  is  more  complicated  than  the  closely 
related  case  of  the  non-resonant  double  Hopf  bifurcation  and  contains  three  different 
types  of  co-dimension  1  singularities  and  4  different  types  of  co-dimension  2 
singularities.  Thus,  with  the  help  of  the  results  presented  in  this  paper,  one  can  apply 
the  analysis  of  van  Gils  et  al.  [11]  to  any  physical  problem  exhibiting  generalized  Hopf 
bifurcation  with  non-semisimple  1:1  resonance.  Furthermore,  it  has  been  shown  by 
Hale  [13]  that,  for  systems  with  linear  operators  whose  superdiagonal  terms  are 
equal  to  1,  an  appropriate  scaling  can  be  used  to  obtain  the  averaged  equations.  In 
the  final  section,  the  averaged  equations  up  to  the  second  order  approximation  are 
obtained  and  compared  with  the  normal  form  equations. 
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II.  Background  and  Notations 

The  problem  of  interest  in  this  paper  is  a  4-dimensional  one.  However,  we  shall 
keep  the  analysis  as  general  as  possible  for  the  time  being.  Consider  a  dynamical 
system  governed  by  autonomous  differential  equations  in  Cn, 

y  =  A  (p)y  +  f(y,  p)  (4) 

where  f:  Cn  -4  Cn  is  a  Cr  vector  field,  r  £  2,  A  is  an  nxn  complex  matrix,  x=0  is  the 
trivial  solution  of  Eq.  (4)  for  all  values  of  p  (i.e.,  flO  ,  p)  =  0)  and  the  nonlinear  vector 
function  can  be  represented  as 


=  f2  (y,p)  +  f3  (y,n)  +-fk  (y,p)  +  •  ■  •  •  (5) 

Here,  we  have  expressed  the  nonlinear  terms  as  a  formal  power  series  of 
homogeneous  terms  with  degree  denoted  by  the  superscripts.  We  define  Hr  to  be  the 
linear  space  of  homogeneous  vector  polynomials  of  degree  k  in  n  variables  with  range 
Cn.  Let  (ej,  ea, ....  ej  denote  the  basis  of  Cn  and  y  =  ,  y2  •  •  yj  be  the 

coordinates  with  respect  to  this  basis.  Thus,  an  element  Pfanjof/f*  can  be 
represented  in  the  form  of  vector-valued  monomials  as 

f*(y,p)  =  X  4GN*K  =  X  4m 0*)^^  ,  |m|  =  k 

t  «,m 


=  X  X  m.  (yr1  y” 2 ” •  ynn)  «. 

*  |  m|  »  k  ' 

with  dim  (ff  (y,p)j  =  (n+k-1)!  /  E(n-1)!  k!]  and  dim  {#£)  =  n  •  dim  |ff  (y,p)|.  Now  that  a 
formal  set-up  for  representing  Eq.  (4)  has  been  obtained,  we  can  consider  the  problem 
of  reducing  Eq.  (4)  to  the  normal  form 

x  =  A(p)x  +  g (x,  p),  g(x,p)  =  g2(x,p)  +  g3(x,p)  +  -.g1‘(x,p)  +  ...  (7) 
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which,  as  stated  previously,  is  in  a  simpler  form  than  Eq.  (4)  and  has  all  the  essential 
features  of  the  flow  near  the  equilibrium  point  of  the  original  system.  The  formal 
solution  of  this  problem  consists  of  determining  near  identity  coordinate 
transformations 

y  =  x  +  h  (x)  ,  h  (x)  =  h2  (x)  +  h3  (x)  +  ■  •  •  hk  (x)  (8) 

where  x  €  fl,  and  Q  is  a  neighborhood  of  the  origin  of  Cn,  such  that  the  analytic 
expressions  of  ffy,p)  are  simplified  to  yield  g(x,p).  Once  again,  f* ,  gkand  hk  are 
homogeneous  vector  polynomials  of  degree  k  and  belong  to  //k  •  Assuming  the  normal 
form  reduction  up  to  order  k-1  has  been  performed,  differentiating  Eq.  (8)  gives 
y  =  [i  +  D*hk  (x)]x 
and  substituting  in  Eq.  (4)  yields 

x  =  [i  +  Dx  hk  (x)]  [a  (x  +  hk  (x))  +  f(x  +  hk  Cx )}] . 

Making  use  of  the  fact  that,  for  x  e  ft, 

[i  +  Dxhk(x>]  1  =  I  -  Dxhk(x)  +  o(|x|2(k*1)) 

results  in 

x  =  Ax  +  fix)  +  f®(x)  +  •  •  •  f^'Hx) 

+  [fk(x)  +  [Ahk(x)-Dxhk(x)Ax]|  +  o(|x|k+1).  (9) 

It  is  worth  noting  that  the  transformation  of  degree  k  does  not  affect  the  normal  form 
of  order  (k-1)  but  does  affect  the  terms  of  order  k  and  higher.  The  task  now  is  to 
select  hk(x)  so  that  the  terms  of  degree  k  in  the  brackets  are  as  simple  as  possible. 
Examining  the  terms  of  degree  k  in  Eq.  (9)  and  comparing  with  those  of  Eq.  (8)  yields 

A  hk(x)  -  Dx hk  (x )  Ax  +  ^(x)  =  g'Hx)  (10) 
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and  fJ (x)  =  g^ (x)  for  j  =  2,3,  •  •  •  k-L  Introducing  a  linear  operator  La  defined  by 

LAhk  =  [hk.  Ax]  =  Ahk(x)  -  Djjh1' (x )  Ax  . 

Eq.  10  can  be  rewritten  as 

-  LAhk(x)  =fk(x)-gk(x)  =  qk  (x ) .  (11) 

The  above  equation  is  called  a  homological  equation.  La  :  -*  is  called  the 

homological  operator  and  is  linear  in  the  space  c,  homogeneous  vector  polynomials  of 
degree  k.  Equation  (11)  is  to  be  solved  for  hk  (x). 

Let  us  denote  Rk  as  the  range  of  La  and  let  be  any  complementary 
subspace  to  Rk  in  ff„  .  H*  can  be  decomposed  as  follows 

ffj  -  ®  wlS  ,  ts2.  (12) 

Thug,  for  each  Ax)  e  there  exists  i)k  (x)  e  and  gk  (x)  €  such  that  any 

given  homogeneous  polynomial  of  degree  k  can  be  written  as 

&(x)  =  g Hx)  -f  q*(x) 

and  the  suitable  transformation  h^x)  is  obtained  from 

-LAhk(x)  =  q^).  (13) 

Since  the  choice  of  complementary  space  is  not  unique,  neither  is  the 
transformation  hk(x)  or  the  normal  form  gk(x).  This  nonuniqueness  was  resolved  by 
Elphick  et  al.  [7]  through  a  particular  choice  of  inner  product.  As  in  [7]  (refer  also  to 
Helgason  [14]),  we  can  introduce  an  inner  product  in  Hn.  To  this  end,  we  introduce  a 
differential  operator  associated  with  an  arbitrary  f*(x)  e  as 
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ffrli-  I  4 


\m 


|m|-k  ^  \d*/ 


e; 


/af  _ 


ami 


aro” 


w  a* 


m. 


dx 


»* 


Then,  for  ff(x) ,  g.k(x)  in  H^,  the  scalar  product  is  given  by 


dx 


(ff(x) ,  g‘(*>)  -4(9)  if  (X)|  t  =  I  I  44-  ~ 

Xm0  lal.k  10- k 


dx 


V 


x»0 


It  is  clear  that  the  only  terms  that  will  survive  are  those  for  which  a  and  p  coincide, 


i.e. 


(f?(*), gfc)}  *  X  m!  =mi!  m2! 

m4c 

Thus,  the  inner  product  in  Hk  is  defined  as 


m. 


*yn  Syn  • 


°  iml  lm|- k 

Using  this  inner  product,  we  can  define  the  anoint  operator  (La)*  as 

f^Ah\x),  i =  (hk  (x) ,  L\f(x))H> 
and  making  use  of  the  fact 

(hk<Ai) ,  f*<4H,  «  (hk&),  fk(A*x))w 
Elphick  et  al.  [7]  has  shown  that 


(14) 


ker(LA-)  =  ker  (La)’.  (15) 

Since  H £  is  a  finite  dimensional  space,  ker  (La)*  is  an  orthogonal  complement  of  Rk 
the  elements  of  which  we  are  free  to  choose.  Equation  (12)  may  then  be  written  as 

Hl  =  ^  ©  ker(LA*)- 


(16) 


Now,  considering  the  linear  equations  in  Hk,  we  have 
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-  LAhk(x)  =  T)k  (x)  ,  La*  gk(x)  -  0 


(17) 


and  the  solvability  condition 

(nk<x)  ,  g k  <*>)».*  *  °*  (18) 

The  normal  form  and  explicit  formulas  for  the  coefficients  can  then  be  calculated 
using  Eqs.  (17)  and  (18).  It  is  important  to  note  that  this  normal  form  depends  on  the 
matrix  A  and  the  choice  of  complementary  space  W^.  Once  the  functions  flKx)  are 
known,  the  above  method  can  be  applied  to  calculate  both  tfKx)  and  gKx).  A  recursive 
algorithm,  similar  to  that  of  Chow  and  Hale  [4],  can  also  be  employed  to  compute  the 
order  nonlinearities  flKx)  given  all  transformations  h(x)  and  normal  forms  g(x)  up 
to  order  k-1.  Both  methods  have  been  employed  independently  herein  to  calculate  the 
normal  form  coefficients  which  are  given  explicitly  in  the  Appendix. 


III.  Normal  Form  for  Non-semisimple  Case 
For  the  non-semisimple  case,  the  normal  form  calculations  are  not  as  easy  as 
in  the  case  of  a  diagonalizable  linear  operator.  However,  the  calculations  can  be 
simplified  using  certain  well  known  results  in  Lie  algebra.  These  will  be  introduced  as 
we  proceed  through  the  calculations  of  the  normal  form  for  the  generalized  Hopf 
bifurcation. 

Given  a  finite  dimensional  vector  space  V  over  the  complex  numbers  C  and  a 
space  L  of  linear  transformations  of  V  onto  itself,  one  can  define  the  Lie  bracket  by 
the  formula 


[P,  Q]  =  (P*  Q-Q*  P)  e  L  for  P,  Q  e  L  . 


(19) 
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Then  L  becomes  a  Lie  algebra  and  we  say  P  commutes  with  Q  iff  [P,Q]  =  0.  The 
result  that  is  of  importance  to  us  is  the  Jordan  decomposition  theorem  which  states 
that  for  any  A  e  L  there  exist  S  and  N  such  that 


A  =  S  +  N  and  [S,N]  =  0  (20) 

where  S  is  semisimple  (diagonalizable)  and  N  is  nilpotent.  Moreover,  these 
decompositions  are  unique  and 


ker  A  =  kerS  kerN.  (21) 

In  the  calculation  of  normal  forms  for  generalized  Hopf  bifurcation  with  non¬ 
semisimple  1:1  resonance,  the  linear  operator  of  interest  takes  the  form 


A  = 


‘  ico  1 
0  id) 

0 

ico  0 
0  ico 

0 

[0  1 
0  0 

0 

0 

-io)  1 

0  -ico  . 

= 

0 

4co  0 

0  -ico  . 

+ 

0 

0  1 

0  OJ 

=  S  +  N 


(22) 


and  [S,N]  =  0.  In  addition,  the  homological  operator  for  any  two  matrices  A  and  B 
also  satisfies  the  relation  [La  ,  Lb]  =  I^A.If.  This  implies  that  the  Lie  brackets  of  Ls, 
Ln  and  Ls*,  Ln*  also  commute.  Thus,  the  ker(LA*),  which  is  needed  for  the 
calculation  of  the  normal  form,  is  given  by 


ker  (La„)  =  ker  Ls*  n  ker  1^* . 

It  is  worth  pointing  out  that  the  above  results  can  also  be  obtained  using  the 
arguments  given  in  Meyer  [15].  Furthermore,  the  normal  form  g(z),  given  in  Eq.  (17), 
commutes  with  elements  of  the  Lie  groups 

G  =  (esA*|se  R)  and  S1  =je,s|s€  r) 
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and  the  normal  form  i9  said  to  have  G  -  equivariance  and  a  simpler  S*  -equivariance, 
respectively.  Since  the  proofs  of  these  results  are  similar,  only  that  of  S1  - 
equivariance,  i.e. 

g(esS£)  =  e^g(%) 

will  be  given  here.  To  this  end,  consider  z  =  e®8  %  and  g (z)  =  g  (e*3  4)  •  Taking  the 
total  differential  of  gfz)  w.r.t.  the  variable  s  yields 

=  D  g  (z)  Sz . 
ds 

Now,  using  the  fact  that  the  normal  form  is  such  that  g  e  ker  (La*)  =  ker(Ls*) n 
ker(LN*)  and  S*  =  -  S,  we  have 

Dx g  (z)  &  -  Sfe(z)  =  0. 

Combining  the  above  two  equations  yields  an  O.D.E  for  g  (z) 

«  Sg (z) ,  seR 
ds 

whose  solution  can  be  written  as 

g  (z)  =  e^g  (z;  8  =  0)  =  e^g  (£)  .  (23) 

This  proves  the  S1  -  equivariance.  The  G  -  equivariance  can  be  proven  similarly  by 
replacing  S  by  A*  in  the  above  steps. 

1.  Linear  Algebraic  Calculation  of  the  Normal  Form  Coefficients 
Now  we  calculate  the  normal  form  and  appropriate  expressions  for  the 
coefficients  of  this  normal  form.  To  this  end,  consider  the  homological  equation 
-  La hk  (x)  =  f*(x). 
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It  is  easy  to  show  that  for  the  semisimple  S  with  eigenvalues 
X, ,  i  =  1,2  ■  •  n ,  LAh^Gc)  reduces  to 

Lshk(x)  =  X  [(m ,  x)  -  Xt]xmef 

•,|m|  ■  k 

and 

(m,  X)  -  X,  =  0;  s  =  l,2,--,n;  |m|  £2 

is  called  the  resonance  condition.  The  ker  (Ls*)  is  determined  by  the  appropriate 
combination  of  m's  which  satisfy  the  above  condition.  The  resonance  condition  for  the 
problem  under  consideration  can  be  expressed  as 

i<o(mi  +  m2  -  m3  -  014  - 1)  =  0 ,  mi+  m2+  003+  m4  =  k  for  a  =  1,2 

and 

•ia>(m3+  m4-m1-m2-l)  =  0,  m1+m2+m3+m4  s  k  for  8  =  3,4. 

Since  mj  £  0  and  integer,  it  is  obvious  that  k  is  always  odd  and  the  above  conditions 
yield 


(mj  +  m2)  =  k*  *  ,  (m3  +  m4)  =  for  s  =  1,2 

i  mt 

and 

K+m^  =  ,  (m3+m4)  =  kil  for  s  =  3,4 . 

Thus,  the  non-zero  nonlinear  normal  form  exists  only  for  k  =  3,5...  .  However,  the 
original  quadratic  nonlinear  terms  can  contribute  to  the  cubic  terms  as  a  result  of  the 
nonlinear  transformation  as  will  be  seen  in  the  subsequent  section.  Calculation  of  the 
coefficients  of  the  leading  order  normal  form  (k=3)  is  of  concern  in  this  paper.  Thus,  in 
an  80  dimensional  basis,  only  24  vectors  lie  in  the  ker  (Ls*)  and  can  be  written  as 

(xf X3) e, , (xi x2 xs) es * U2X3) e»  >  (x?x4)eg,  (xix2x4)e,,  (xlxje,  fors  =  l,2 
(x3 i) e» » (*3 x4 x  1) e, ,  (x4x J e, ,  (x^x2)e#,  (x3x4X2)e,,  (xjxje,  for  8  =  3,4. 
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The  action  of  Ln*  on  these  bases  can  be  represented  by  a  24x24  matrix  of  the 


form 

010110' 

0  0  2  0  1  0 
0  0  0  0  0  0 
0  0  0  0  1  0 
0  0  0  0  0  2 
0  0  0  0  0  0- 

where  I  and  0  are  6x6  identity  find  zero  matrices.  The  8-dimensional  null  space  of  the 
above  matrix  can  be  easily  computed.  Making  use  of  this,  the  basis  of  ker  (La*)  can 
be  written  as 


CO  0  0 

-I  C  0  0 

0  0  CO 
0  0  -I  C 


where  C  = 


X1(X1X4-X2X3) 

X2(xix4-xax3) 

0 


\ 

2 

*1  *3 

Xi2X3 

1  ° 

Xi2X4 

1 

*1*2X3 

1  *12*3 

r 

0 

9 

0 

’  I  0 

1 

0 

0 

1  0 

0  1 

,  0  , 
*3  (*2*3  -  X1X4) 

1  »  ■ 

0 

0 

Xq^Xi 

0 

0 

Xq  X, 

1  »  ! 

OOO 

*4  (*2*3- X1X4)  1 

*3  1 

2 

3  1 

*3  *2  ' 

*3*4*1 

i  *3  *1  / 

It  is  worth  noting  that  the  first  4  basis  vectors  are  complex  conjugates  of  the  last  4, 
as  expected.  Since  any  linear  combination  of  these  vectors  spans  the  null  space,  we 
can  manipulate  the  given  basis  such  that  the  resulting  normal  form  is  as  simple  as 
possible.  This  manipulation  is  performed  as  follows:  the  second  basis  element  is 
replaced  by  the  vector  obtained  by  subtracting  the  third  basis  element  from  the 
second  and  the  sixth  basis  element  is  replaced  by  the  vector  obtained  by  subtracting 
the  seventh  basis  element  from  the  sixth.  This  procedure  yields  the  new  2nd  and  6th 
bases  as 
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0  1 

0  i 

*i  <*1*4- *»*»>!  and 

0  \ 

0  ( 

0  1 

0 

X3  (XgXg-XjX^J 

Thus,  the  normal  form  for  the  generalized  Hopf  bifurcation  with  1:1  resonance  can  be 
written  as 

(ijH’okJ  (zj)+  ,a,(z‘Zl)*  a2*Zl?2 '  ZlZ2)l  (zj)+  NwHbsl** '  ®  ) 

(24) 

where  aj  *  Cj  +  idj  ,  bj  *  ej+ifj  ,  j  =  1,2.  In  the  above  equation,  we  have 
replaced  (xi,  X2,  X3,  X4)  by  (zi,  Z2,  Zi,  z2).  Thus,  the  second  and  third  equations  can 
be  obtained  by  conjugating  the  above  equations. 

While  calculating  the  coefficients,  we  shall  assume  that  the  original  system 
contains  both  quadratic  and  cubic  nonlinearities.  Thus,  for  the  problem  under 
consideration  in  this  paper 

f*(y)  =  X  X  ^i/n2^i3ni4  (yTV? 2y3  3y?4)  .  dim  (fff)  *  4C  (25a) 

s*l  |ml  ■  2 

f3(y)  =  X  X  Cia»2W4Wly2ay33y44)e.*  di“(^4)  *  80 *  (25b) 

|m|»3 

We  have  shown  in  the  previous  section  that  ker  (La*)  =  10)  for  k  =  2  since 
ker(Ls*)={0}.  Thus,  all  the  quadratic  terms  given  by  expression  (25a)  can  be 
eliminated  and  the  transformation  which  performs  this  reduction,  obtained  by  the 
matrix  representation  of  ker(LA*),  is  given  by 
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where 


r _ 

hl;m 

f|m 

r  b  i 

0  O' 

0  B  0  0 

h2;m 

f|m 

0  0  B  I 

hsjn 

L  0  C 
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and  hj'm  and  ff^n  are  vectors  of  dimension  10.  Since  B  is  nonsingular,  it  is  easy  to 
calculate 


him  *  B1^  ,  him  =  B'^ffpn-hiJ 

and  h3;m  and  h4?n  are  the  conjugates  of  hi^  and  h2?n ,  respectively.  The 
complete  expressions  for  h|m  and  hijm ,  are  given  explicitly  in  the  Appendix.  As 
these  transformations  annihilate  all  of  the  quadratic  nonlinearities  in  the  given 
system,  they  alter  the  terms  of  order  3  and  above.  We  denote  the  new  coefficients  of 
the  cubic  nonlinearities  as 

=  C +  C  C-HCCCC) 

where  f|m  are  the  original  coefficients  of  the  cubic  nonlinearities,  and  f«-m  are  the 
coefficients  of  the  new  cubic  terms  generated  while  eliminating  the  original  quadratic 
nonlinearities.  The  coefficients  are  indeed  functions  of  the  coefficients  of  the  original 
quadratic  nonlinearities  as  one  would  expect.  Now,  the  normal  form  for  the  leading 
nonlinearity  is  given  by  Eq.  (23)  and  is  defined  in  the  space  complementary  to  R4.. 


1  6 


The  coefficients  ai ,  a2 ,  bi ,  b2  and  their  conjugates  are  calculated  using  the 
solvability  condition  of  Eq.  (18).  The  first  4  coefficients  are 

ai  »  ^  (3  fi;2010  +  f|;1110  +  fi^OOl}  +  ai 

a2  =  ^  (  2f^;2001  -  2f2;0210  -  $;1110  +  f|;110l)  +  a2 

bi  =  ft;20io  +  hi 

b2  =  4  (fl^OlO  *  f|;1110  +  3  f|,200l)  +  hz 

4 

where  the  expressions  for  ai,  S2,  bi  and  b2  in  terms  of  the  coefficients  of  the 
quadratic  nonlinearities  are  given  in  the  Appendix.  The  remaining  four  coefficients 
are  obtained  by  conjugation  of  the  above  expressions,  i.e., 
a3  =  ai,a4  =  a2,h3  =  bi  and  b<  =  l>2- 

2.  Recursive  Calculation  of  Normal  Form  Coefficients 
This  approach  is  based  on  a  series  of  papers  by  Ponce,  Gamero  and  Freire  [8, 
9, 10, 17]  which  are,  in  turn,  implementations  of  a  method  of  Chow  and  Hale  [4,  Chap. 
12]  which  employs  a  technique  developed  by  Deprit  [6]  using  Lie  transforms  to 
determine  the  normal  form. 

In  order  to  remain  consistent  with  the  literature,  the  following  notation  will  be 
used:  define  F*,  tA,  G*  e  Hn  by 

Fk  =  (k  - 1)!  f*(y) ,  (/k  =  (k- 1)!  hk(x) ,  Gk  =  (k-1)!  g\x) . 

The  first  step  is  a  rescaling  to  isolate  the  homogeneous  terms  of  degree  k.  Letting 
x— >ex  and  y— >ey  for  e  e  R,  the  original  system  of  Eq.  (4)  becomes 
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(26) 


[k‘fl 

W  M 
W  4  M 

(fJ)  ^  M 

(fJ)  F*  F*  F?  [f£] 

(!)  !  :  !  -  ( :]. 
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The  terms  in  round  brackets  are  from  the  original  system  and  those  in  square 
brackets  are  the  final  normal  form.  Each  term  in  the  triangle  depends  on  those 
immediately  to  the  left  and  above.  The  indexing  scheme  used  here  is  different  from 
that  in  the  above  references;  the  superscript  k  refers  to  both  the  order  of  the 
monomials  in  the  vector  and  the  row  in  which  it  appears  in  the  Lie  triangle  and  the 
subscript  refers  to  the  column  in  which  it  appears  in  the  Lie  triangle. 

The  recursion  operates  across  rows  of  the  Lie  triangle  from  left  to  right.  As  an 
example  of  what  occurs  during  a  recursion,  consider  the  fifth  row.  F1  is  a  vector 
containing  the  order  5  terms  in  the  original  system.  To  generate  F2,  Fl  is  added  to  the 
sum  of  the  terms  in  column  1  above  F3  combined  with  the  appropriate  U  i's; 
F4  x  U2,  x  U3,  f\  x  U4,  Fj  x  U 6.  To  generate  f|,  F2  is  added  to  the  sum  of  the 
terms  in  column  2  above  F2  combined  with  the  appropriate  U  i's; 
F2  x  U 2,  F%  x  U3,  F%  x  U4.  This  process  is  continued  until  F3  is  reached  at  which 
time  the  normal  form  has  been  obtained.  What  is  happening  as  the  recursion  moves 
across  the  Lie  triangle  is  the  accumulation  of  the  order  5  contributions  of  the  near 
identity  transformations  of  orders  2  to  5.  In  column  2,  the  contributions  from 
substituting  the  transformations  into  the  original  equations  are  collected.  In 
succeeding  columns,  the  contributions  from  the  interaction  of  new  terms  generated 
by  the  transformation  and  the  subsequent  transformations  are  collected  until  finally 

[k-M 

in  column  5  one  has  the  order  5  terms  of  the  normal  form.  The  coefficient  \  j  -  2  /  which 
appears  in  the  sum  is  a  counting  term  analogous  to  the  binomial  coefficient  in  the 
binomial  theorem. 

Now  rewrite  Eq.  (11)  using  the  new  notation 
-L AUk=Fk-Gk 

where  F*  is  a  vector  of  the  order  k  monomials  resulting  from  the  near  identity 
transformations  up  to  order  k-1.  If  Eq.  (11)  is  rewritten  as 


then 


Fk  =  Gk-LAtfk 


Pr<4«rLA/k  =  Gk 


prqjgkF  =-L AU 
However,  if  (11)  is  written  as 


Gk  =  Fk+LAl/k 


and  it  is  noted  that  F^  x  t/k  =  (Ax)  x  t/k  =  LAt/k,  then 
Gk»Fk  +  Fj  x  f/k. 

Now  consider  the  recursion  (29).  The  only  time  U  k  will  appear  is  when  Z  =  2,  in 
which  case  (29)  can  be  written 


f]  x  C/k 


f£=f2  +  f}  x  C/k . 
For  Z  =  3,  Eq.  (29)  can  be  written 


Fj  =  F2k  +  Fjxl/k+  ilf-l)****** 


Fk*F2 +FjXt/k. 


So,  for  any  Z,  2  :£  Z  :£  k, 


Fk  =  F k  +  Fj1  x  (/k 


It  is  easy  to  see  that  the  Fk  obey  the  recursion  relations 
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This  recursion  is  identical  to  Eq.  (29)  except  for  F2  ,  there  the  F^  x  U*  term  was  left 

k 

out  (and  thus  will  not  appear  in  any  of  the  subsequent  Fl ).  Thus 
Fkk  =  Fk  =  Gk-LA(/k. 

So  Eq.  (30)  can  be  used  to  determine  the  order  k  normal  form,  and  if  Eq.  (31)  is 
written 

J/k  =  -L>oiR;Fk 

the  order  k  near  identity  transformation  can  be  obtained. 

In  order  to  continue  on  to  higher  order  terms  in  the  normal  form,  it  is  necessary 
to  convert  the  F^  s  into  if*  s.  This  is  accomplished  using  the  following  correction 

Fk  =  Fk  -  projR* Fk ,  /  «  2, . ..,  k. 

IV.  Dominant  Normal  Form 

In  order  to  study  perturbations  of  a  vector  field  with  linear  part  given  by  the 
non-semisimple  matrix  A,  we  consider  the  universal  unfolding  of  the  linear  vector  field 
Ax  used  in  van  Gils  et  al.  [11] 

A(X)  =  f1+a  .  1  ),  ^(a,^,^),  ae  R,  pspi  +  ipje  C 

14  »  +  «/  (32) 
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This  may  be  calculated  explicitly  using  the  homological  equation  (10)  applied  to  first 
degree  polynomials  hl(x),  in  the  same  manner  as  (23)  for  cubic  polynomials.  The 
unfolding  parameters  X.  are  found  in  terms  of  the  original  linear  coefficients  to  be 


a“  2 {^i;10o°  +  4)010°)  ^  = 


2  “zyiwj  -  r  2^000 

The  above  unfolding  of  A(X)  may  also  be  found  from  the  viewpoint  of  versal 
deformations  of  matrices,  as  in  Arnold  [1],  allowing  for  rescaling  of  time. 

Now,  making  the  observation  that  zi  =  0  implies  Z2  =  0  and  the  normal  form 
commutes  with  S,  we  choose  a  transformation  as  in  van  Gils  et  al.  [11] 
zj  =  re®,  Z2  =  re*®w,w  =  u  +  iv,  =  cot  +  9 
which  yields  three  real  equations  independent  of  the  phase  variable  9 

r  =  r  [a  +  u  +  r2  (ci  +  2d2  v)] 
u  =  pi  -  u2  +  v2  +  r2  (ei  +  2f2  v) 
v  =  p2  -  2uv  +  r2  (fi  •  2e2  v) 
and  0  =  v+  r2  (d1-2c2v). 

In  order  to  "blow  up"  the  dominant  terms,  we  rescale  the  above  variables  as 
r  =  e r ,  u  =  e u  ,  v  =  £  v,  et  =  t,  a  =  e  a ,  pi  =  e2  pi ,  42  =  e2  P2-  Introducing  r  2  =  p 
and  dropping  the  hats,  we  have,  in  new  time, 


and 

where 


p'  =  2p  (a  +  u)  +  e2ci  p2  +  0  (e2) 

u'  *  pi  -  u2  +  v2  +  eip  +  e2f2  pv  +  0  (e2) 

v’  =  p2  -  2uv  +  fip  -  e2e>2  pv  +  O  (e2) 

0’  =  v+  edi  p  +  0(e2) 

ci  =  Re  (ai)  ,  ei  =  Re  (bi) ,  e2  =  Re  fo) 
dx  =  ImiaJ  ,  fx  =  Im(bJ,  =  Im^). 


(33) 
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V.  Averaged  Equations 


In  this  section,  we  shall  demonstrate  the  relationship  between  second  order 
averaging  and  normal  forms  for  the  nilpotent  case  under  consideration.  To  this  end, 
we  make  use  of  the  scaling  suggested  by  Hale  [13]  for  linear  operators  whose 
superdiagonais  are  equal  to  1.  In  order  to  make  the  calculations  less  cumbersome,  we 
only  consider  cubic  nonlinearities  and  the  nonlinear  system  can  be  written  as 


y  =  A(p)y  +  F°(yi,  y3)+  F1(y)+  F2(y)+  F*(y2,y4)  (34) 

where  A  is  as  given  in  Eq.  (22)  and  the  nonlinearities  of  degree  3  can  be  written  in 
terms  of  the  original  notation  as 

=  X  j4ooo  y? +  *82010  y?y3 +  4<>20  y^ +  *!po3o  y^)  e8 

8*1 

pl  =  X  (**2100  y?  y2  +  **2001  yfo  +  4no  y^  +  4on  y^ 

•-1 

+  *8)0120  y2y3+  4>021  y^y4)e. 

f2  =  X  (4200  y^+  4  002  +  4ioi  y^*  4m  y2y3y4 

+  *8  0210  *8)0012  ytfl) e. 


8*1 


4  . 

^  *  X  (*$0300  & +  42oi  &4  +  4io2  yA  +  4kx»  y4| e*  • 

8*1 

In  order  to  bring  the  above  equations  into  "standard  form”,  we  make  use  of  the  scaling 
suggested  by  Hale  [13],  which  is  in  line  with  that  of  van  Gils  et  al.  [11], 
yi  =  oci ,  y2  =  £2x2 ,  Y3  =  oc3 ;  y4  =  &A 
ana  transform  Eq.  (34)  to  new  variables  z  by  means  of  the  transformation 
~  rr  »  Ai*2  -  V 


x,  =  zje*"'  ,  xW2  =  ,  j  =1,2. 


This  procedure  yields  a  set  of  equations  in  standard  form  to  O(e^)  as 


z  a  eX°(z,z,t)  +  e2X1(z,z,t)  ,  za(zi,z2) 


(35) 


where 


z2 

-rl  _ 

e  -totFi(zi,zi,t) 

.e-**FS(»i,*i,t). 

f  — 

.  e  -*ntF2(z,z,t)  . 

(36) 


and  the  z  equations  are  obtained  by  conjugating  Eq.  (35).  Now,  applying  the 
averaging  procedure  up  to  the  second  order  yields 


z  =  eM  (X°(z,z,t)!  +  e2M 


ax 


ax°  ,-tt  aw  ^  aw 


w  +  w  - 

dz  dz  dz 


_  X  +  X1  (z,z,t)j 
dz  | 

(37) 


where  M  is  the  averaging  operator  defined  as 

rT 


M  (•  )  =  lim  if  (•  )dt 
t  t-*-  T  J0 


and 


wia 


-r 


X°dt+  c(z,z)  ,  X  (z,z,t)  a  X°(z,z,t)  -  M  (X°(z,z,t)} 


i.  e. 

a  Cj(z,z)  and  W^(z,z,t)  ak(z,z,t)  +  c^z) 

with  k(z,  z,  t)  defined  as 


k(z,Z,t)  a  ?2$(M  Z1  e2k#t"  2i^  *2*020  Z1Z1  e  2kDt*  *20030  Z1  e  4* 


4kot 


where  c  is  an  arbitrary  vector  function  of  z  and  z.  The  choice  of  c  is  made  such  that 
the  normal  form  coincides  with  the  resulting  second  order  averaged  equations.  We 
have  made  two  observations  concerning  the  product  terms  within  the  second  curly 
bracket  of  Eq.  (37).  Note,  in  Eq.  (36),  that  X2  is  only  a  function  of  z\  and  z\  and 
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Xi  =  0  ,  Xj  sb  0  Thus,  the  second  order  contribution  from  k(z,z,t)  is  identically  zero. 
The  second  order  contributions  to  the  averaged  equations  are 

e2M  (c2-L(c1)+x}(z,z,t)j 
e2M  (2^10 

C1Z1Z1  +  4;2010  ci  Zj  -  L  (C2H  X2  (z,z,t)| 


z,  Z? 


,  r  M  «.)_  .  ac.)  #s  2 ».)y  at.)  1® 

where  3^  **  +  ^  W*i  2i  ’ 

Comparing  terms  of  like  order  in  the  averaged  and  normal  form  equations,  the 
appropriate  choice  of  the  vector  c  is  given  by 


Cj(z,z)  s=  ax  Zg  ,  C2(Z,Z)  =  0t2  zf  Zj  . 

Equating  coefficients  yields 

al  4^010"  **2  =4(^010  *  4^001  "  4;lllo]  *  al  4^010 '“2  ”4  [*12010  "  4#00 1  ‘  4a  110 


It  is  obvious  that  otj  must  be  real.  Choosing  aj  to  be  identically  zero,  Oj.  is  obtained 
as 


Thus,  the  averaged  equations  are 

i1  *  e  Zj  +  Aj  (zj  Zj)  Zj  +  CKe3) 

®  £bj  (zjZj)  Z1  +  e2  [1)2  (z,%)  +  -  h,)  (z^)]  z1  +  CKe3) .  (38) 

The  second  pair  of  equations  are  obtained  by  conjugating  Eq.  (38).  As  before,  we 
introduce  the  universal  unfolding  defined  by  matrix  A(X)  (see  Eq.  (32))  into  Eq.  (38), 
use  the  transformation 
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(39) 


zi  a  re'8  ,  Z2  =  reiew,  w  =  u  +  iv 
and  rescale  the  variables  as 

et=t,  \i i  =  z2VLXt  U?  =  e2 £2 -  a  =  m.  (40) 

After  substituting  Eqs.  (39)  and  (40)  into  Eq.  (38)  and  dropping  the  hats,  we  have  the 
averaged  equations  in  terms  of  p  =  r2  as  expressed  in  Eq.  (33).  Thus,  one  can 
conclude  that  the  dominant  terms  of  the  scaled  normal  form  equations  (33)  agree 
completely  with  those  of  the  averaged  equations. 
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Appendix 

2 

The  transformations  h  ,  i  =  1,2,  which  eliminate  the  quadratic  terms  are: 

ipn 


h2  S-L(  ?  +i(1)f2  ) 

1,2000  ^  \  2;2000  1;2000| 


**1:1100  *  ^3  [®(*  ^1;2000  +  *2;110o)  +  1  (^-^OOO  +  ^^l.lioo)] 

^1;0200  “  ^4  ['  ®2;2000  +  0)2  (*  *l;1100  +  *2;020o)  +io)(*  2*i;2000  +  2*2;1100  +  °>2^1;02CK>)  ] 


1:1010  “  ^ 
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ABSTRACT 

A  generalized  four  dimensional,  nonlinear  and 
non-autonomous  system  is  studied.  The  effect  of 
periodic  parametric  excitations  is  examined  on  systems 
that  exhibit  Hopf  bifurcation  with  one  to  one  resonance 
along  with  subharmonic  external  resonance.  The  linear 
operator  is  assumed  to  have  a  generic  nonsemisimple 
structure.  In  this  case,  the  dimensionality  of  the  system 
can  net  be  reduced  despite  the  presence  of  an  S1 
symmetry.  However,  the  system  is  simplified 
considerably  by  reducing  it  to  the  corresponding  four 
dimensional  normal  form  equations.  The  local  behavior 
of  the  equilibrium  solutions  is  studied  along  with  their 
stability  properties.  Several  codimension  1,  2  and  3 
bifurcation  varieties  are  observed.  Some  of  the  global 
bifurcations  that  are  present,  can  be  associated  with  the 
Bogdanov  Takens  and  10,  +i,  -i)  bifurcation  varieties. 
The  numerical  results,  obtained  by  using  AUTO  and 
CHAOS,  indicate  the  existence  of  homoclinic  orbits 
along  with  the  period  doubling  behavior  which  leads  to 
chaos. 

*•  INTRODUCTION 

The  objective  of  this  study  is  to  investigate  the  effect 
°f  periodic  parametric  excitations  on  systems  that 
•xhibit  Hopf  bifurcation  with  1:1  resonance.  For  this 
Purpose  we  consider  a  generalized  four  dimensional, 
nonlinear,  non-autonomous  system.  It  is  assumed  that, 
■t  some  critical  parameter  values,  the  linear  operator 
c°ntains  two  coincidental  purely  imaginary  eigenvalues 
*hich  generically  lead  to  a  non-semisimple  structure, 
nder  these  conditions,  the  system  defined  on  a  four 
'mensiona!  center  manifold  is  described  by  a 


5-parameter  family  of  normal  form  equations. 

One  to  one  resonant  Hopf  bifurcation  has  been 
studied  by  several  authors  in  the  past  (Caprino  et  &1. 
(1984);  Meyer  (1984);  Krupa  (1986))  and  recently  by  van 
Gils  et  al.  (1990),  who  considered  3-parameter  unfolding 
of  the  vector  field  singularity.  The  work  of  van  Gils  et  al. 
considered  a  system  similar  to  one  presented  in  thiB 
paper  but  without  parametric  excitation,  and  the 
associated  normal  form  was  reduced  to  three 
dimensional  first  order  equations.  Here,  several 
codimension  1  and  codimension  2  bifurcation  varieties, 
in  addition  to  some  interesting  periodic  behavior,  were 
examined.  Furthermore,  it  was  shown  that  the  normal 
form  of  the  fourth  order  autonomous  system  has  an  S1 
symmetry,  and  hence  it  is  possible  to  reduce  the 
dimension  by  one,  bringing  the  system  to  a  set  of  three 
equations  in  an  appropriate  co-ordinate  frame.  In  this 
paper,  we  consider  the  corresponding  nonautonomous 
generic  system.  As  in  van  Gils  et  al.  the  structure  of  the 
equations  in  the  present  case,  can  be  simplified 
considerably  by  bringing  them  into  their  normal  form  or 
by  using  the  method  of  averaging.  Here,  a  modified 
version  of  the  normal  form  of  Namachchivaya  et  al. 
(1992)  is  considered  which  takes  into  consideration  the 
subharmonic  resonance  with  respect  to  the  forcing 
frequency,  in  addition  to  the  1:1  internal  resonance. 
However  the  disadvantage  in  such  a  system  with 
parametric  harmonic  excitations  is  that  the 
dimensionality  of  the  system  can  not  be  reduced. 
Moreover,  in  addition  to  the  traditional  3  unfolding 


parameters  for  the  autonomous  system,  one  has  to 
introduce  2  more  parameters,  namely  the  amplitude  of 
the  forcing  and  a  detuning  parameter  which  represents 
the  deviation  of  the  excitation  frequency  from  twice  the 
natural  frequency.  One  expects  very  rich  local  dynamics 
in  such  a  system  due  to  the  presence  of  5  parameters. 
For  the  problem  under  consideration,  there  is  an  S1 
symmetry  and  the  normal  form  equations  can  be  further 
simplified.  The  equilibrium  points  of  this  simplified 
fourth  order  system  extend  to  the  periodic  orbits  that  are 
the  relative  equilibria  of  the  original  four  dimensional 
normal  form. 

The  system  of  equations  discussed  here  arises  in  a 
variety  of  physical  problems  which  have  1:1  internal 
resonance  under  parametric  periodic  excitation.  The 
problems  of  parametric  excitation  of  nonlinear 
dynamical  systems  are  of  importance  in  several 
branches  of  engineering  such  as  vibrations  of  beam 
structures  under  dynamic  loads,  flow  induced  vibrations 
and  control  systems.  The  major  part  of  this  analysis 
deals  with  the  local  behavior  of  the  equilibrium  solutions 
and  their  stability  properties.  We  also  study  the  global 
bifurcation  set  associated  with  Bogdanov-Takens  and 
(0,  +i,  -i)  bifurcation  varieties.  In  addition,  we  observe 
period  doubling  bifurcation  which  leads  to  chaos.  Due  to 
the  complicated  nature  of  the  computations  involved,  we 
have  used  symbolic  computations  extensively.  The 
analytical  solutions  have  been  compared  with  those 
obtained  numerically  using  AUTO  [Doedel  (1986)},  a 
general  purpose  software  package  for  the  bifurcation 
analysis  of  differential  equations  and  the  maps.  We 
have  also  used  CHAOS  [Aronson  (1991)},  which  is  a 
SUN-based  software  program  for  numerically 
simulating  the  nonlinear  systems.  The  analysis  of  the 
complete  global  behavior  is  in  progress  and  we 
anticipate  presentation  of  these  results  in  the  near 
future  [Namachchivaya  and  Malhotra  (1992)]. 


2.  STATEMENT  OF  THE  PROBLEM 

We  consider  the  following  system  on  a  four 
dimensional  center  manifold: 


quadratic  and  cubic  nonlinearities,  and  the  linear 
operator  A  =  Dxf(0)  has  two  equal  purely  imaginary 
nonzero  eigenvalues  along  with  their  conjugates. 
Furthermore,  the  double  eigenvalues  are  assumed  to  be 
non-semisimple,  and  the  time  dependency  on  f(x,t)  is 
given  explicitly  in  terms  of  parametric  harmonic 
excitation.  Under  these  conditions,  the  system  (2.1)  can 
be  expressed  as: 

y  =  Ay  +  Pf  B  y  Cos  twf  t)  +  f  (  y,y  )  (22) 

where  y  e  C4,  yt,3  =  xj  ±  i  x^,  ya,4  =  X3  ±  i  X4  i  gf  and  a*  are 
the  amplitude  and  the  frequency  of  the  parametric 
excitation  and  B  is  a  constant  4x4  matrix.  f(y,y)  is  the 
nonlinear  term.  The  linear  operator,  A,  has  the 
following  non-semisimple  form: 

iw  1  00] 

*  0  i  to  0  0 

0  0  -iw  1 
0  0  0  -i  w 

where  i  0)  is  the  pure  imaginary  eigenvalue.  In  the 
subsequent  sections  we  study  the  stability  and  the 
bifurcation  behavior  of  the  system  described  by  Eq  (  2.2). 


3.  TRANSFORMATION  TO  STANDARD  FORM 

In  order  to  study  the  dynamic  behavior,  it  is 
desirable  to  reduce  the  system  under  consideration  (2.2) 
to  its  normal  form.  The  normal  form  of  the 
nonautonomous  (periodic)  linear  part  of  Eq.  (2.2)  is 
obtained  via  a  periodic  transformation  y  =  v  +  H(t)  v  (See 
Namachchivaya  and  Malhotra  (1992)),  as: 


v  =  Av  +  ^ 
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1  U),  t 

b13Vle 

1  (Or  t  _  1  t 

b23  Ve  +  b13v2  e 


b3i  vi  e 


i  0%  t 


+  h.o.t. 


•  1  <J)t  t  -  1  Wf  t 

b41  v,  e  +  b3,  v2  e 


x  =  f(x,t)  (2.1) 

where  x  €  R4,  is  the  state  vector  and  f:R4  — *  R4  ,  is  the 
non-autonomous  vector  field  which  is  analytic  in  its 
arguments.  We  assume  that  the  system  (2.1)  has  both 


We  define  a  new  time  i  such  that  t  =  uif  t  ,  where 
wf  =  wo  (  l  -  e  X  ).  and  X  is  the  detuning  parameter.  In 
addition  to  the  internal  1:1  resonance,  we  consider  the 
case  of  subharmonic  parametric  resonance  (k  =  u  wo 
=  1/2).  Without  loss  of  generality,  letting  cuo  =1.  and  using 
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i  i,o  new  time  derivative  with  respect  to  t,  Eq.  (3.1)  takes 
the  following  form: 


I  V,  Mr  .  _ 

v,  =  — +  v.,  +  ■  b.,  v,  e  i  + 

i  2  2  2  13  1  i 


f  \±1 1  +  X  v2  +  X  bJ3  Vj  eM 


+  h.o.t. 


(3.2) 


v2  =  i  ^2  +  7*  ( b23  V1  e' 1  +  b13  v2  e' 1 ) 


*[“-  +  ^<*b23vi  e,T  +  Xbl3v2  e"! 


+  h.o.t. 


The  universal  unfolding  of  A,  obtained  from  a  four 
parameter  versal  deformation  given  by  Arnold  (1983),  i8 
expressed  as: 


A  (a) 


(iu>+a)  1 

p  (ico+a) 


(3.5) 


where  o  =  (a,  p)  and  p  =  (pj  +  i  P2).  van  Gils  et  al.  (1990) 
interpret  a  as  a  real  crossing  parameter  and  p  is 
interpreted  as  a  complex  splitting  parameter  for  the 
eigenvalues  of  the  linear  operator  of  A.  In  order  to  get 
rid  of  the  time  dependent  terms  in  Eqs  (3.3),  we  make  the 
following  transformation. 


The  prime  denotes  differentiation  with  respect  to  i,  and 
the  other  two  equations  can  be  written  as  complex 
conjugates  of  these  equations. 

The  next  step  is  to  include  the  normal  form  of  the 
nonlinear  terms  in  Eqs  (3.2).  For  a  detailed  discussion  of 
the  calculation  of  the  nonlinear  normal  form  with  a  four 
dimensional  non-semisimple  linear  operator,  one  is 
referred  to  Namachchivaya  et  al.  (1992).  Following  this 
procedure  yields  the  complete  normal  form  of  the 
original  system  of  Eq.  (2.2)  and  is  expressed  as  following: 


,(iV2) 


V3,<  =  U3.«  e 


C-1V2) 


(3.6) 


The  nonlinear  terms  remain  invariant  under  this 
transformation  due  to  S1  symmetry  and,  in  order  to 
distinguish  the  dominant  terms  in  Eqs  (3.3),  we 
introduce  the  following  scaling: 


ui,3  =  e  2i,3  > 


n  o 

“2,4  =  £  224-  M-f  =  £  h  ,  and 
a  =  £  P,  p  =  t2  u 


(3.7) 


vj  =  [^-  +  v2  +  |f  bi3vi  e'T+aj  v)2v1  +  a2v12v2 

-  a2  V1  v2  V1  1  +  £  X  l  2  +  v2  +  2  bl3  V1  e 
+  a,  Vj2  Vj  +  a2  Vj2  v2  -  a2  Vj  v2  Vj  }  +  h.o.t. 

(3.3) 
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+  a2  v,  v2  v2  -  a2  Vj  v22  +  b2  Vj2  v2  ]  +  h.o.t. 


Now  consider  the  universal  unfolding  of  the  linear 
operator  A  as  given  by  (2.3).  The  two  dimensional 
analog  of  A  can  be  written  as: 


A  = 


ICO  1 

0  iw 


(3.4) 


On  applying  the  universal  unfolding  given  by  (3.5),  the 
transformation  given  by  (3.6)  and  the  scaling  given  by 
(3.7),  the  normal  form  Eqs  (3.3)  become 
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PZj  +  Zj  + 


iX  Zj 


E2[gbj3  z,+  Xz2  +  aj  Zj2  z,]  +  o(e3) 
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Z,  =  E 
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where  aj,  bj,  82>  b2  can  be  expressed  in  terms  of  the 
original  nonlinear  coefficients  of  Eq.  (2.2).  It  is  worth 
mentioning  that  these  equations  up  to  0(ej)  terms,  have 
also  been  obtained  by  using  the  method  of  averaging,  as 
shown  in  Namachchivaya  et  al.  (1992).  The  above 
mentioned  reduced  system  has  5  parameters.  The  three 
unfolding  parameters  (P,  ui.u?)  control  the  behavior  of 
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the  eigenvalues  of  the  linear  operator  A,  while  g  and  X  redundant  stability  conditions  are  expressed  as: 
control  the  extent  of  the  forcing  in  terms  of  its  amplitude 

and  the  deviation  from  the  subharmonic  excitation  Tl  :  P<  0 

frequency.  With  the  help  of  these  parameters,  it  is 

possible  to  explore  the  local  dynamics  of  the  system  in  T2  (u,  +  +  f-*  +  >  gb"  ,4  3, 

the  neighboring  regions.  In  order  to  explore  the 

dominant  dynamics,  the  time  t  is  rescaled  as  t  =  t/e  so  T3  :  Uj  <  '  U30  ui  *  uw> 

that  in  the  slow  time  Eqs  (3.8)  take  the  following  form: 

where. 


*1  *  p  zt  +  z2  +  +  0(e) 

*• 

(3.9) 

Xj  *  u  +  P  z,  +  gb  Zj  +  +  bj  Xj2  Zj  +  0(e) 

where  gb  *  g  b;o  ,  and  the  dot  denotes  differentiation  with 
respect  to  the  slow  time. 

4.  STABILITY  OF  THE  TRIVIAL  SOLUTION 

It  is  dear  that  z  =0,  z  =  0  is  the  trivial  solution  of  the 
normal  form  equations,  and  as  a  first  step  we  consider 
their  stability.  Transforming  Eqs  (3.9)  to  real  variables 
(*1,  yi,  X2,  y2>  by  means  of  the  usual  transformations 


U‘°  *  2  (v  *  2pi'  =  -  1^.  vi3(j  = 
and 

4(?2  (4(3*  +  ft2  X2  +  gb2)  (44) 

40  '  (4P2  +  X2) 

Thus  T2  and  T3  completely  determine  the  stability  of  the 
linear  operator  L  once  Tl  is  ensured.  The  typical 
stability  region  for  the  trivial  solution  of  (4.2)  is  shown  in 
Fig.  (1)  in  (U|,  uj)  parameter  space,  where  T2  represents 
the  rerion  outside  the  circle  centered  at  (-U10,  -U;o)  with 
radius  gb.  while  T3  represents  the  region  on  the  left  side 
of  the  parabola  described  by  u;2  +  U30  t>i  =  uio 


Xj  «  (zj  +  zt)/2 ,  yl-  (zj  -  zt)  /(2i)  and 

Xj,  *  (zj  +  Zj)/2  ,  y2=  (z2-z2)/(2i),  5.  EXISTENCE  OF  NON  TRIVIAL  SOLUTION 


yields  the  following  set  of  equations: 

0 
0 

c2k3  +  xiy2)  -  d2(*2y1+yi3) 
c2(*2yi+  y3) +  ^l*3  +  *iy2) 

where 


In  order  to  examine  the  existence  of  the  non-tnvia! 
equilibrium  solutions,  we  need  to  transform  Eqs  3  9  into 
a  convenient  coordinate  system.  The  structure  of  these 
equations  suggest  the  coordinate  transformation  of 
(zi,  Z2)  to  the  new  coordinates  (r,  0,  u,  v)  by  using 

z\  =  re  Z2  -  re  (u  +  iv)  (5  j) 

where  (r,  u,  v)  e  R3,  and  0  is  the  phase  variable  This 
leads  to  the  following  set  of  four  coupled  equations. 


p 

-A/2 

1 

0 

X/2 

P 

0 

1 

K  +  gfe) 

*  U2 

p 

-A/2 

K ' 8b) 

X/2 

P 

p  =  2p(p  +  u)  +  0(E) 
u  =  Uj  +  (v2  •  u2}  +  c2  p  +  gb  cos20  +  0(e) 
v  =  u2  -  2  u  v  +  d2  p  -  gb  sin20  +  0(e) 


(5.2) 


and  \n  and  t)2  are  the  real  and  imaginary  parts  of  the 
scaled  complex  splitting  parameter  u,  and  C2  and  d2  are 
the  real  and  imaginary  parts  of  the  nonlinear  coefficient 
bi  in  Eqs  3.9.  For  the  linear  operator  L,  the  non 


9  =  v  +  +  0(e) 

2 


where  p  =  r2  >  0.  This  coordinate  transformation 

involves  a  singular  change  of  the  coordinates,  which 


where 


rest  nets  the  appiicability  of  Eqs  (5.2)  near  p  =  0,  but  this 
possibility  has  already  been  considered  in  the  previous 
section.  In  the  absence  of  parametric  excitation  (i.e  . 
gb  =  0  and  A  =0).  Eqs  *5.2)  reduce  to  a  set  of  three  coupled 
equations  (in  p,  u,  v),  independent  of  the  phase  variable  0 
and  a  fourth  equation  describing  the  evolution  of  0. 
which  is  precisely  the  form  studied  by  van  Gils  et  al. 
(1990).  The  parametric  excitation  acts  to  couple  the 
phase  (B)  equation  with  the  other  three  equations.  The 
equilibrium  points  of  Eqs  (5.2)  extend  to  the  periodic 
orbits  of  the  normal  form  of  the  original  system.  Thus 
the  study  of  the  non-trivial  equilibrium  solutions  of  the 
above  system  is  vital  to  our  analysis.  The  dominant  part 
of  Eqs  (5.2)  is 

p  =  2p(p  +  u) 


a  =  (c22  +  d22) 

b  *  2[c2(u,  +  ujo)  ♦  d2  (\>2  ♦  t>2o)3  (5.7) 

c  *  (iij  ♦  \>io?  ♦  ('U2  +  U20)2  *  gb2 

and  uio  and  020  are  defined  in  Eqs  (4.4).  The  non-trivial 
equilibrium  solutions  (5.5)  are  admissible  only  if  p  is  real 
ar  ’  positive. 

According  to  Implicit  Function  Theorem,  these 
equilibrium  solutions  (5.5)  extend  smoothly  to  the 
equilibrium  solutions  of  (5.2)  for  small  e  >  0.  if  the 
Jacobian  matrix  (Jo)  of  (5.5)  is  nonsingular.  The 
Jacobian  matrix  (Jq)  is  given  by: 


li  =  Uj  +  (v2  -  u2)  +  c2  p  +  gb  cos20 

(5.3) 

v  =  u2  -  2  u  v  4  d2  p  -  gb  sin20 

9  =  v  +  ^ 

2 

The  equilibrium  solutions  of  Eqs  (5.3)  can  be  obtained  by 
solving  the  following  algebraic  equations 

P  4  u  =  0 


0 


0 


where 


2p  0  0 

2P  -A  e1 

X  2  p  e2 
0  1  0 


(5.8) 


e,  = +  +  e2  =  2[u, +Ui0  +  C2p) 


Det  [Jo]  =4p  b  4  a  pj , 


(5.9) 


Uj  +  (v2  -  u2)  +  c2  p  +  gb  cos20  =  0 
v2  -  2  u  v  +  d2  p  -  gb  sin20  =  0 


a  and  b  are  given  by  Eqs  (5.7).  We  observe  that  if  the 
(5.4)  discriminant  (D)  of  the  quadratic  equation  (Eq.  (5.6))  is 
nonzero,  then  the  Jacobian  Jo  (5.8)  is  nonsingular,  i.e., 


v  4  A/2  =  0 


D=(b2-4ac)£0  =»  DetIJO]  *  0  except  where  D=0 


Eqs  (5.4)  yield  the  following  equilibrium  points 

-b±Vb2-4  ac 
P  =  - 2 

u  =  -P 

v  =  -A/2 

0  =  1  sin'1 
2 

p  satisfies  the  following  quadratic  equation: 
a  p2  4  bp  4  c  =  0 


9 

u  -  2  u  v  4  d2  p 
gb 


(5.5) 


(5.6) 


This  implies  that  the  equilibrium  solutions  (5.5)  extend 
locally  in  e,  to  the  equilibria  of  (5.2),  except  where  D  is 
zero.  Throughout  this  investigation  the  nondegeneracy 
condition  i.e.,  (C22  +  d22)  >  0  is  assumed  to  be  satisfied. 
Eq.  (5.6)  suggests  that  p  may  have  2,  1  or  0  equilibrium 
solutions  that  satisfy  the  requirements  on  p.  Depending 
upon  the  values  of  the  parameters  we  have  the  following 
scenarios. 

Case  A: 

The  following  should  be  satisfied  for  two 
equilibrium  solutions  of  p  to  exist,  which  would  lead  to 
rpal  and  positive  values  of  p. 

A  l:  Cj  (v,  +  o10)  4  d,,^  4  Ujq)  <  0 
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A.2:  (\>!  +  v10f  +  (u2  +  u2 -  gb  2  >  0 

A.3:  {<^2  (Uj  +  Ujo)-  C2(u2  +  u2(j)  |  <  I  V[c22  +  djj2)  .  gb| 

Case  B: 

If  a  >  0,  and  c  <  0,  Eq.  (5.6)  will  have  only  1  root  that 
satisfies  the  requirement  on  p.  In  this  case,  the  sign  of 
the  coefficient  b  does  not  affect  the  number  of 
equilibrium  solutions  of  p,  but  it  determines  the 
magnitude  of  the  positive  root  of  p.  Thus,  the  following 
condition  needs  to  be  satisfied  for  only  one  root  of  p  to 
exist. 

B-l:  (ux  +  vj  +  (u2  +  Ujq)2  -  gb2<  0 

Elsewhere  in  the  parameter  space,  p  does  not  have 
a  permissible  solution.  Figure  (2)  shows  the  typical 
regions  where  2,  1  or  0  non-trivial  equilibrium  solutions 
of  p  exist  for  different  values  of  the  nonlinear 
coefficients.  For  both  cases  (Fig.  2a  and  b),  there  is  one 
permissible  root  inside  the  circle,  two  such  roots  outside 
the  circle  between  two  tangent  lines  and  no  real  and 
positive  roots  elsewhere.  Furthermore,  the  second 
stability  condition  for  the  trivial  solution  T2,  and  the 
existence  condition  A.2  and  B.l  are  given  by  the  same 
circle  centered  at  (-uio,  -U20)-  For  the  nondegenerate 
case  (a  >  0),  various  possibilities  for  the  number  of  non¬ 
trivial  equilibrium  solutions  are  shown  in  Fig.  (3), 
where  ftp)  is  plotted  vs.  o. 

6.  STABILITY  OF  THE  NON-TRIVIAL  SOLUTION 

In  this  section  we  discuss  the  stability  of  the  non¬ 
trivial  equilibrium  solutions  obtained  in  the  previous 
section.  The  Jacobian  matrix  (5.8)  has  a  dependence  on 
the  equilibrium  values  of  p,  which  in  turn  depend  on  the 
system  parameters  in  a  complex  way  as  given  by  Eqs 
(5.5).  Using  the  Routh-Hurwitz  criteria,  we  arrive  at  the 
following  stability  conditions  that  must  be  satisfied  for 
the  non-trivial  equilibrium  solution  (if  it  exists): 

Nl:  p  <  0 

N2:  (uI  +  u10)  <  2 jp2  +  L.J  .  2c,p 

N3:  {2P(Uj  +  ul0)  +  X  (\>2  u20)  +  2  (2  pt^  +  X  d2)  p)  >  0 

N4:  (2  +  ap)  51  0  (6.1) 


N5:  [80  p6  +  24  p4  X2  +  p2  X*  -  96  p*  u4  -  8  p2  X?  u, 

+  16p2  ut2  -  32  p3  X  u2  -  4  X2u22  -  128  p4  c2  p 

2  2  2  2 

-  32  P  X  Cj  p  -  64  p  u2  d2  p  -  16  X  u2  d2  p 

-  64  p2  d22  p2  - 16  X2  d22  p2  ]  >  0 


where  p  is  given  by  (5.5).  Due  to  the  complex  nature  of 
these  conditions,  it  is  difficult  to  comment  on  the  nature 
of  the  stability  of  the  non-trivial  equilibrium  solution. 
Figure  (4)  indicates  the  several  solution  branches  of  p. 
which  are  given  by: 


P 


Jj_  Vb2  4  ac 
2a  2  a 


(6.2) 


p*  «  1(.  i  +  Vb^ac  )lb<0c>0)  (63; 

p+1  =  |(-  £  +  — 2~a4-  )lb  <  0,c  <  0)  (6.4) 

p  *2  =  1  (-  ^  )  I  b  >  0,c  <  01  (6.5) 


Clearly,  the  lower  branch  (p  )  violates  N4,  which 
requires: 

Thus  p'  is  an  unstable  branch  of  the  non-trivial 
equilibrium  solution.  It  can  be  shown  that  p ♦  and  p*1 
branches  are  unstable  since  they  violate  the  sufficiency 
condition  (N5).  Numerically  it  can  easily  be  checked  that 
part  of  p+2  equilibrium  branch  satisfies  all  five 
conditions  of  stability  (Nl  to  N5).  The  complexity  of  the 
computations  hinders  the  calculation  of  the  stable  part  of 
p+2  equilibrium  branch,  analytically.  Sometimes  the 
fourth  condition  (N4)  can  provide  important  information 
about  the  stability  of  the  non-trivial  equilibrium  solution 
branches. 


7.  BIFURCATION  ANALYSIS 

In  this  section  we  describe  the  bifurcation  behavior 
associated  with  Eqs  (4.1).  The  steady  state  solutions  of 
(4.1)  represent  periodic  solutions  of  the  original  system. 
In  this  system,  several  codimension  1  and  codimension 
2  bifurcation  varieties  are  detected  along  with  a 
codimension  3  singularity.  One  is  referred  to 
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Sumachchivaya  and  Malhotra  (1992)  for  a  detailed 
description  of  the  necessary  conditions  for  the  different 
bifurcation  varieties  to  exist.  Following  is  a  brief 
description  of  each  of  these  bifurcation  varieties  that  are 
detected  in  our  system  (4.1): 

This  bifurcation  variety  occurs  when  one  of  the 
eigenvalues  of  the  linear  operator  L  passes  through  zero. 
When  this  happens,  the  condition  T2  is  violated.  Under 
these  conditions,  the  trivial  solution  loses  its  stability 
and  a  non-trivial  solution  with  p  >  0  bifurcates  from  the 
trivial  solution.  For  fixed  values  of  (J,  X  and  gb,  condition 
T2  can  be  depicted  as  a  circle  (&)  in  the  parameter  space 
(ui  ,U2),  as  shown  in  Fig.  (1).  Along  any  generic  path 
transverse  to  £>,  the  stable  trivial  solution  would  lose  its 
stability  (provided  T1  and  T3  are  satisfied)  through 
simple  bifurcation,  when  the  following  condition  is 
satisfied: 

(u,  +  u10)2  +  (u2  +  u20)2  =  gb2  (7.1) 

When  the  linear  operator  L  has  a  zero  eigenvalue,  the 
system  (4.1)  can  reduced  to  its  corresponding  one 
dimensional  center  manifold.  If  we  fix  (5,  X,  gb  and  Uj  at 
-0.5,  0.1,  0.5  and  0.0,  respectively,  the  one  dimensional 
manifold,  in  terms  of  the  nonlinear  coefficients  (C2  and 
d2)  and  the  perturbation  parameter  t)2  =  v>2  -  v>2cr,  is  given 
as: 


provided  T1  and  T2  are  satisfied.  In  this  case,  the  linear 
operator  L  has  a  pair  of  pure  imaginary  eigenvalues 
with  nonzero  imaginary  part  at  the  critical  value  of  the 
parameters.  This  bifurcation  variety  occurs  along  a 
generic  path  transverse  to  the  curve  1),  which  is  depicted 
as  parabola  in  in  the  parameter  space  (uj ,  \>2),  as  shown 
in  Fig.  (1),  provided  other  stability  conditions  are 
satisfied.  The  dynamics  of  the  system  (4.1)  can  be 
studied  on  a  two  dimensional  canter  manifold.  Once 
again,  for  the  same  values  of  0,  X,  gb  and  V},  the 
corresponding  two  dimensional  center  manifold  can  be 
expressed  as: 

A.  at  \>2cr  =  - 0.705354, 

f  =  r  [(0.5ti1-0.7Ti2)  +  (3.8c2-7.18d2)r2  +  ...] 

6  =  [o.426  +  (-  0.5  Tij  -  0.8  tj2  )-*-(-  4.2  c,,  -  8.45  dj.)  r  2  +  ...] 

B.  at\>2cr~  0.705354, 

r  =  r  [(0.5  rij  +  0.7  r\2 )  +  (2.86  c2  +  4.8  dj)  r  2  +  ...] 

8  =  [o.568  +  (-  0.5  q,  +  0.6  q2  )  +  (-  2.7  4  +  4.1  dj)  r2  +  ...] 

where  Tii  and  02  are  the  perturbations  in  the  critical 
values  of  ui  and  u2.  One  can  easily  observe  that  the 
periodic  equilibrium  branches,  as  described  by  these 
center  manifold  equations,  have  the  same  dynamical 
behavior  as  observed  numerically,  in  the  neighborhood 
of  the  critical  value  of  the  bifurcation  parameter  \>2 
(Fig.  7). 


at  u2cr  =  -  0.484446, 

w  =  v.  ‘2.13  o2-  19.98  n22  +  (c2+  1.76 d2)w2] 

and, 

at  u2cr=  0.484446, 

w  =  w  -  1.49  r\2 ' 6-97  q22  +  (c2  *  18  d2)  w2] 

It  can  be  easily  verified  that  the  nature  of  the  non-trivial 
equilibrium  branches,  as  captured  by  these  center 
manifold  equations,  agrees  qualitatively  with  the 
numerical  results  (Fig.  7). 

7.2.  Hfipf  Bifurcation: 

This  bifurcation  variety  occurs  when  condition  T3  is 
violated,  i.e.: 


7.3.  Bogdanov  Taheaa  Bifurcation: 

Now  we  study  the  possibility  of  the  intersections  of 
two  codimension  1  bifurcation  varieties,  the  simple 
bifurcation  ft  and  the  Hopf  bifurcation  J).  At  this 
intersection  (C2)  we  have  two  eigenvalues  of  the  linear 
operator  L,  becoming  zero,  where  C2  is  defined  by: 

4  &t>,  +  2  X  u2  -  p(*P*  +  X2}  =  0 

(v,  +  vl0f  +  (t>2  +  =  gb2 

These  equations  lead  to  the  following  two  critical  points 
in  (uj  ,t>2)  parameter  space. 


(7.3a) 

(7.3b) 


%>U1  +  U40 


(7.2) 
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t>,CT  =  •  ± 


gb  X 


«a  =-V20± 


10  I - 2 -  2  ’ 

V  (4  P  +  X  ) 

2  gb  P 

V(4P%X2) 


(7.3c) 


At  each  of  these  critical  points  the  linear  operator  can  be 
brought  to  the  following  form: 


of  aa  and  b3  can  be  scaled  to  unity  by  reseating  other 
variables.  For  positive  values  of  U2,  as  \>2  is  varied  from 
zero  for  a  fixed  value  of  uj.  we  obtain  the  phase  portraits 
6,  7,  8,  9,  10  and  11  as  indicated  in  Fig.  (5),  these  show 
various  local  and  global  bifurcations  and  completely 
agree  with  the  numerical  results  presented  in  the  next 
section. 


0  10  0 
0  0  0  0 

0  0  X3  0 

0  0  0  x4 


(7.3d) 


where  X3  and  X4  represent  the  two  stable  eigenvalues  of 
L  at  the  critical  points  (7.3c).  Dynamics  near  this 
nonhyperbolic  fixed  point  can  be  studied  by  examining 
the  equations  defined  on  a  two  dimensional  center 
manifold.  The  procedure  to  reduce  the  four  dimensional 
system  to  the  corresponding  two  dimensional  center 
manifold  is  fairly  systematic.  Once  the  reduction  to  the 
two  dimensional  center  manifold  is  achieved,  we  obtain 
the  corresponding  truncated  normal  form  (Gamero  et 
al.  (1991))  Now  we  compare  our  system  with  a  similar 
system  studied  by  Takens  (1974),  which  is  as  following: 

W1  =  W2 

Wj  =  04  Wj  +  ctj  w2  +  a3  w,3  +  b3  W]2  w2 

where  a*  and  a2  are  the  unfolding  parameters,  wj  and 
W2  are  the  two  dimensional  center  manifold,  and  33  and 
b3  are  the  nonlinear  coefficients.  Considerable  research 
has  already  been  done  on  this  bifurcation  variety. 
Guckenheimer  and  Holmes  (1983)  discuss  the  details  on 
the  global  bifurcation  behavior  associated  with  this  class 
of  bifurcation. 

If  we  fix  p,  X  and  gb  at  -  0.5, 0.1  and  0.5  respectively, 
the  linear  operator  (4.3)  has  double  zero  eigenvalues  at 
the  following  critical  parameter  values  of  u,  an  Uj  : 

A.  vt<r  =  0.29725,  u*,  =  0.44752 

If  q,  and  qj  are  the  small  perturbations  in  and 
Ujor .  then  the  corresponding  values  of  ct| ,  a2 ,  aa  and  b5 
are  given  in  terms  of  qi  and  q2  and  the  original 
nonlinear  coefficients  C2  and  d2,  in  the  following 
manner: 

dj  =  -  0.11  q,- 1.09  q2,  0^  =0.97  qt  - 1.4  q2 

a3=-  0.134  %  - 1.34  ,  b3  =  6.26  Cj  ♦  8.28 

For  C2  =  1-  and  d2  =  1„  83  <  0  and  l>3  >  0,  the  magnitude 


B.  Vl<r  =  0.19775,  =  -  0.54752 

The  corresponding  values  of  ai  ,  a; ,  a,  and  bs  are 
given  as  following: 

Oj  =  0.1  qt  +  0.9  q2  ,  a2  =  0.8  q,  -  1.  q2 
a3  =  0.1  c2  +  0.7  dj  ,  b3  =  2  6  c2  -  4.15  d2 
For  C2  =  1.  and  d2  =  1.,  33  >  0  and  b3  <  0.  the  magnitude 
of  a3  and  b3  can  be  scaled  to  unity  as  before  by  rescaling 
other  variables.  For  the  negative  values  of  02,  once  again 
as  l>2  is  varied  from  zero  for  a  fixed  value  of  uj ,  we  obtain 
the  phase  portraits  5,  4,  3,  2  and  1  as  shown  in  Fig.  (5j. 
Thus  combining  these  two  cases  we  can  can  construct 
partially  the  bifurcation  diagram  for  the  original 
problem. 

Depending  upon  the  values  of  a3  and  b3.  two  distinct 
cases  of  unfolding  (i.e.,  I:a3>0,  b3  <  0  and  II:  83  <  0, 
b3  <  0)  are  possible  for  this  bifurcation  variety.  As 
mentioned  before,  these  two  cases  are  completely  studied 
by  Takens  (1974)  and  Guckenheimer  and  Holmes  (1983). 
The  other  two  cases  (Ic:  83  >  0,  (53  >  0  and  llc:  a^  <  0, 
1>3  >  0)  can  be  constructed  from  I  and  II  by  reversing  the 
sign  of  W2,  a2  and  time.  Table  1  shows  the  possible  cases 
of  unfoldings  that  can  occur  for  all  possible 
combinations  of  C2  and  d2,  near  the  two  codimension  2 
bifurcation  points  A  and  B. 


Table  1 


Case 

Possible  unfolding 
(A) 

Possible  unfolding 
(B) 

C2  >  0,  d2  >  0 

II' 

I.  lc 

C2  =  0,  d2  >  0 

II' 

I 

C2  >  0,  d2  <  0 

I,  Ic,  IK 

Ic.  II' 

C2  =  0,  d2  <  0 

I 

II' 

C2  <  0,  d2  >  0 

I,  II.  IIC 

I.  II 

C2>0,  d2  =  0 

11' 

I' 

C2  <  0,  d2  <  0 

I 

II,  11' 

C2  <  0,  d2  =  0 

I 

II 
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7.4.  Simple  and  Hopf  Bifurcation  Variety: 

In  this  case,  the  linear  operator  has  a  simple  zero 
eigenvalue  along  with  a  pair  of  pure  imaginary 
eigenvalues  for  some  critical  value  of  the  parameters. 
Following  are  the  requirements  for  the  simultaneous 
existence  of  simple  and  Hopf  bifurcation  varieties: 

(i>i  +  u10f  +  (u2  +  u20f  =  gb2  (7.4a) 

4  (3  -  2  X  v2  -  p  (20  p2  +  X2)  =  0  (7.4b) 


Uj  <  (3  P2  +  — 


These  conditions  lead  to  the  following  critical  point  in 
(ui,t>2)  parameter  space: 


222 
P  (4p  +  X  ) 


»  5P  •  V  +  P*. 


^2cr  =  *  P^  +  2P 


2  2  2 
P  (4  p  +  X  ) 


At  the  critical  point,  the  system  can  be  reduced  to  its 
three  dimensional  center-manifold  and  the  linear 
operator  assumes  the  following  form: 


0  -co  0 
-  <0  0  0 

.000 


The  structure  of  the  linear  operator  completely 
determines  the  normal  form  for  the  linear  and 
nonlinear  parts  of  the  three  dimensional  center- 
manifold  equations.  If  we  fix  p,  X  and  gb  at  -0.1, 0.1  and 
0.5  respectively  then  the  resulting  normal  form,  in 
cylindrical  co-ordinates  is  given  as: 


r  =  r  (  Hj  +  au  r2  +  a,2  z2  )  +  h.o.t. 
z  =  z  (  p2  +  a21  r2  +  z2  )  +  h.o.t. 


0 

where, 


to  +  h.o.t. 


Pi  a  0.04  q,  +  0.95  q2 
P2  *  2.2  q,  -  5.6  q2 
aM  =  0.58  Cj  +  3.42  d2 
a12  =  -  5.64  c2  +  68.15  d2 


a2i  a  4.72  C2  -  42.66(12 
=  55.84  4  -  142.25  cL, 
to  =  0.64 

pi  and  P2  are  the  unfolding  parameters  that  are 
expressed  in  terms  of  q*  and  q2  which  represent  the 
perturbations  in  the  values  of  Uj  and  \>2  from  their 
critical  values,  i.e., 

qi  =  ui  -  uicr,  qa  =  \>2  -  U2CT 
where  t>icr  =  -  0.175,  \>2cr  =  0.455,  and  ajj's,  the 
nonlinear  coefficients  of  the  the  normal  form,  are  given 
in  terms  of  nonlinear  coefficients  cj  and  d2- 

We  notice  that  the  6  equation  is  decoupled  from  the 
r  and  z  equations,  so  we  can  study  the  planar  system 
(r,z)  independent  of  0  and  later  interpret  the  results  for 
the  full  three  dimensional  system  (r,0,z).  A  rescaling 
procedure  can  be  applied  to  the  first  two  equations  of  7.4f 
to  bring  these  to  the  form  described  in  Guckenbeimer 
and  Holmes  (1983).  After  rescaling,  these  equations  up 
to  third  order  are  expressed  as: 

r  =  r(p,+  r2  +  b  z2 ) 

(7.4g) 

zsztpj  +  cr^  +  dz2) 
where  d  =  ±  1,  b  =  and  c  *  -^21_ 
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This  system  has  been  studied  by  Takens(1974),  Langford 
and  Iooss  (1980)  and  Holmes  (1980).  In  this  work,  we  use 
the  same  numbering  scheme  as  given  in  Guckenheimer 
and  Holmes  (1983).  They  classify  12  different  cases  of 
distinct  unfoldings.  Table  2  shows  the  possible  cases  of 
unfoldings  that  can  occur  for  all  possible  combinations 
of  C2  and  d2,  the  nonlinear  coefficients  of  our  original 
normal  form  equations. 

Table  2 


Case 

Possible  unfoldings 

C2  >  0,  d2  >  0 
C2  =  0,  d2  >  0 

la.  Ib.  II,  III.  V,  Via,  VIb.  Vila,  VHb 

II.  Via,  VIb 

C2  >  0,  d2  <  0 

III.  Vila,  V 

C2  a  0,  d2  <  0 

III,  Vila,  Vllb 

C2  <  0,  d2  >  0 

II,  Via,  VIb 

C2  >  0,  d2  a  0 

III.  Vila.  Vllb 

C2  <  0,  d2  <  0 

II.  III.  IVa,  IVb.VIa.VIb.VIIa.VIIb.VIII 

C2  <  0,  d2  *  0 

II.  Via.  VIb 

If  we  start  with  system  (7.4g)  and  consider  the  case  where 


C2  *  da  =  1-0,  one  can  easily  notice  that  the  resulting 
system  corresponds  to  Via  class  of  unfolding.  The 
bifurcation  set  and  the  associated  phase  portraits  for  this 
class  of  unfolding  are  discussed  in  the  literature.  When 
viewed  for  full  three  dimensional  flow  of  Eq.  (7.4f),  these 
results  qualitatively  match  with  the  phase  portraits 
obtained  by  numerically  integrating  the  original  normal 
form  equations  (Eqs  (4.1))  near  the  relevant  parameter 
values. 


7JS.  Double  Hoof  Bifurcation  Variety: 

Now  we  consider  the  case  where  the  linear  operator 
L  has  two  pairs  of  identical  purely  imaginary 
eigenvalues.  The  conditions  for  this  codimension  2 
bifurcation  variety  can  be  simplified  to  the  following; 

X  =  0,  P  =  0,  vi  <  0,  U2  =  ±  gb  and  Ui2  +  V22  >  gb2  (7.5a) 


Under  these  conditions,  the  linear  operator  takes  the 
following  form: 


0  0 

0  0 

(u,  +  gb)  -  v2 
u2  (u,  -  gb) 


1 

0 

0 

0 


0 

1 

0 

0 


(7.5b) 


Whether  Lh  will  have  the  semisimple  structure  or  non- 
semisimple,  depends  upon  the  value  of  gb.  For  the 
stability  of  the  trivial  solution,  we  require  P  <  0,  if  P  =  -  e, 
where  0  <  e  «  1,  the  typical  bifurcation  diagram  for  this 
case  is  presented  in  Fig.  (6).  Along  a  path,  transverse  to 
®2.  the  trivial  solution  loses  its  stability  through  double 
Hopf  bifurcation  as  P  -»  0.  In  fact  t>2  =  ±  gb  is  the 
degenerate  form  of  the  1)  curve  as  P  ->  0  and  X  -»  0. 

I ):  v2  ■  '  uaoyi  +  v40  (7.5c) 

As  p  -*  0  and  X  -»  0,  we  have  ujo  0  and  U40  -*  gb2, 
thus  U2  *  ±  gb. 

7  A  Bogdanov  Tokens  and  Hopf  Bifurcation: 

Along  with  the  assortment  of  codimension  1  and 
codimension  2  bifurcation  varieties,  we  may  also  have 
the  case  where  the  linear  operator  has  a  pair  of  zero 
eigenvalues  along  with  a  pair  of  pure  imaginary 
eigenvalues.  This  codimension  3  bifurcation  variety  has 
not  been  studied  in  detail.  The  conditions  for  the  linear 
operator  L  to  go  through  this  bifurcation  variety  are  as 
following: 


P  <  0,  X  *  0,  ui  =  p2  and  \>2  =  ±  gb  i7.6) 

This  is  also  a  degenerate  bifurcation  variety  (£3).  Along 
a  path  through  £3,  the  trivial  solution  may  lose  its 
stability  through  (0,  0,  ±  i)  bifurcation  variety. 

The  trivial  solution  may  lose  its  stability  in  one  of 
the  many  ways  mentioned  in  this  section,  giving  rise  to  a 
non-trivial  steady  state  branch  or  a  periodic  equilibrium 
branch  of  solutions.  The  non  trivial  equilibrium  may 
also  lose  its  stability  through  a  secondary  bifurcation. 
The  secondary  hopf  bifurcation  occurs  when  condition 
N5  is  violated.  In  this  case  a,  secondary  equilibrium 
branch  of  periodic  solutions  is  generated.  These  periodic 
branches  may  lose  their  stability  after  going  through  a 
period  doubling  or  homoclinic  bifurcation,  which  may 
lead  to  chaotic  behavior.  The  other  bifurcation  varieties 
(of  higher  codimension)  are  also  possible,  but  those  are 
not  considered  here. 


8.  NUMERICAL  SIMULATIONS 

Numerical  computations  were  performed  in  order 
to  verify  the  results  obtained  in  the  previous  sections  and 
to  study  in  more  detail  the  various  bifurcation  varieties 
and  the  structure  of  the  periodic  orbits.  The  numerical 
simulations  were  carried  out  using  AUTO,  a  software 
package  for  continuation  and  bifurcation  problems  in 
ordinary  differential  equations.  The  structure  of  the 
periodic  orbits  was  studied  in  detail  using  CHAOS,  a 
versatile  software  for  simulating  nonlinear  systems. 
These  numerical  simulations  were  very  helpful  in 
confirming  the  bifurcation  behavior  obtained  through 
analysis  and  in  providing  a  substantial  clue  for  the  new 
phenomena,  which  were  global  in  nature  and  not 
detected  through  the  theoretical  analysis.  The 
parameter  regions  where  period  doubling  or  chaos  may 
occur,  were  located  as  a  result  of  these  simulations. 

Figure  (7)  shows  the  1-parameter  bifurcation 
diagram  for  the  normal  form  Eqs  (4.2).  In  this  diagram 
the  coordinate  yi  is  shown  vs.  u a ,  as  u,  is  varied  between 
-  1.0  and  1.0,  while  P,  v,,  X  and  gb  are  fixed  at  -  0  5,  0.0. 
0.1  and  0.5,  respectively.  The  nonlinear  coefficients  c., 
and  da  are  fixed  at  1.0.  By  changing  the  sign  of  either  c< 
or  dj  ,  the  orientation  and  the  nature  of  the  stability  of  the 
non-trivial  equilibrium  and  periodic  branches  change, 
but  the  bifurcation  behavior  remains  qualitatively 
similar  for  codimension  1  bifurcation  variety.  In  these 
diagrams  the  solid  (dashed)  curves  represent  the  s  able 
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lun-.iabl*')  stationary  solutions,  solid  (open)  circles 
represent  stable  (unstable)  periodic  orbits  and  solid 
(open)  squares  represent  Hopf  (Simple)  bifurcation 
varieties.  As  \>2  is  increased  from  -1.0,  the  unstable 
trivial  equilibrium  branch  first  gains  stability  through  a 
Hopf  bifurcation  (at  label  2),  which  results  in  a 
equilibrium  branch  of  stable  periodic  solutions.  The 
trivial  solution  loses  its  stability  through  a  simple 
bifurcation  and  results  in  an  unstable  non-trivial  steady 
state  equilibrium  branch  (at  label  3).  As  ua  is  increased 
further,  the  trivial  branch  of  equilibrium  goes  through  a 
Simple  and  a  Hopf  bifurcation  and  results  in  a  stable 
steady  state  branch  (at  label  4)  and  an  unstable  periodic 
branch  of  equilibrium  solutions  (at  label  5).  This  non¬ 
trivial  equilibrium  branch  loses  its  stability  through  a 
secondary  Hopf  bifurcation  and  gives  rise  to  a  secondary 
branch  of  stable  periodic  solutions.  This  matches 
identically  with  what  has  already  been  shown  in  Fig.  1, 
2a  and  4  of  the  previous  section. 

The  periodic  branch  (through  label  2)  approaches 
the  unstable  non-trivial  equilibrium  branch  and  loses  its 
stability  by  going  through  a  homoclinic  bifurcation.  The 
transition  from  the  stable  periodic  orbit  to  the  homoclinic 
orbit  is  clearly  visible  in  Fig.  (8  a  &  b).  The  labels  in  this 
diagram  (Fig.  8)  correspond  to  the  labels  in  Fig.  (7).  The 
unstable  periodic  branch  (at  label  5)  rises  up  and  tends  to 
approach  the  unstable  non-trivial  equilibrium  branch, 
but  it  merges  with  the  secondary  branch  of  stable 
periodic  solutions  (through  label  11).  A  very  interesting 
phenomenon  is  observed  just  before  the  merger.  The 
unstable  periodic  branch  (through  label  5)  goes  through 
a  cascade  of  period  doubling  bifurcation,  which 
eventually  leads  to  chaos  (Fig.  (9)).  Precisely  at  the 
merging  point  a  homoclinic  orbit  is  observed  (Fig.  (8d)). 
This  homoctinic  bifurcation  is  due  to  the  proximity  of 
Bogdanov-Takens  bifurcation  point.  The  lower  branch  of 
the  non-trivial  equilibrium  solution  (through  label  4) 
also  loses  its  stability  through  a  secondary  bifurcation, 
which  results  in  a  symmetry  breaking  stable  periodic 
branch,  which  rises  up  and  goes  through  a  period 
doubling  bifurcation  near  the  unstable  trivia!  solution 
branch  and  ends  with  the  period  becoming  very  large. 

To  study  the  period  doubling  sequence  and  the 
possibility  of  chaos,  computations  were  performed  using 
CHAOS  near  the  point  where  the  change  in  stability 
occurs  for  the  global  periodic  branch  (from  label  5). 
These  computations  were  performed  for  the  same  values 
of  parameters,  except  v,  was  set  at  0.1  to  get  closer  to  the 


region  where  the  Bogdanov  Takens  bifurcation  occurs. 
At  u2  =  0.3589,  one  stable  periodic  orbit  is  seen  ,  which 
goes  through  a  period  doubling  bifurcation  at  u2  =  0.35882 
(Fig.  (9a  &  b)).  These  periodic  branches  go  through  a 
cascade  of  period  doubling  bifurcation  sequences,  and 
finally  become  chaotic  at  Uj  =  0.35875  (Fig.  9e  &  f).  This 
chaotic  behavior  quickly  changes  into  a  regular  two 
period  orbit  structure,  which  also  goes  through  period 
doubling  sequence  and  leads  to  another  chaotic  attractor 
(at  Uj  =0.35855). 

Now,  in  order  to  study  the  bifurcation  behavior  near 
(0,  +i,  -i)  singularity,  we  fix  P,  vi,,  X  and  gb  at  -  0.1,  -  0.3, 
0.1  and  0.5  respectively,  while  fixing  C2  and  d2  at  1.0.  As 
noted  in  the  previous  section,  this  singularity  occurs  at 
(t),„  =  -  0.175,  v2„  =  0.455).  In  this  case  (Fig.  (10)),  the 
steady  state  bifurcation  behavior  is  similar  to  what  had 
been  observed  earlier,  but  the  nature  of  the  periodic 
solutions  which  result  due  to  secondary  hopf  bifurcation, 
is  different.  Both  the  periodic  branches  (through  labels  4 
and  6)  lose  their  stability  after  going  through  a  period 
doubling  bifurcation.  Figure  (lla-e)  shows  the  phase 
portraits  near  the  va  values  where  the  period  doubling 
sequence  occurs,  and  finally  the  flow  becomes  chaotic  at 
u,  =  -0.352. 

To  conclude,  the  numerical  results  confirm  the  rich 
variety  of  local  bifurcation  behavior  observed 
analytically.  Codimension  1  and  codimension  2 
bifurcation  varieties  are  observed  in  numerical 
simulations.  In  addition  the  numerical  simulation 
exhibits  breaking  of  the  homoclinic  orbits,  change  in 
stabilities,  rich  periodic  behavior  and  the  period 
doubling  cascade  leading  to  chaos.  These  bifurcation 
diagrams  give  a  clue  to  study  the  global  bifurcation 
behavior  that  is  associated  with  the  Bogdanov-Takens 
and  (0,  +i,  -i)  bifurcation  varieties,  which  were  described 
in  the  previous  section. 

9.  CONCLUSIONS  AND  DIRECTIONS  FOR  FUTURE 

RESEARCH 

In  this  study,  we  analyze  the  stability  and 
bifurcation  behavior  of  a  parametrically  excited,  four 
dimensional  nonlinear  system  under  the  combined 
conditions  of  one  to  one  internal  and  subharmonic 
parametric  resonance.  The  stability  properties  of  the 
trivial  and  non-trivial  solutions  of  the  5-parameter 
family  of  normal  form  equations  were  investigated  and 
various  codimension  1,  2  and  3  bifurcation  varieties  were 
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located.  The  reduced  fourth  order  system  may  have  2,  1 
or  0  non-trivial  equilibrium  solutions  depending  upon 
the  parameter  values.  These  steady  state  solutions 
extend  to  the  periodic  orbits  which  are  the  relative 
equilibria  of  the  original  four  dimensional  normal  form. 
The  global  bifurcation  behavior  associated  with  double 
zero  and  (0,  +i,  -i}  eigenvalues  is  also  investigated  and  it 
matches  identically  with  the  numerical  results.  The 
complete  analysis  of  the  global  bifurcation  behavior  for 
this  system  is  currently  in  progress.  The  numerical 
results  indicate  homoclinic  bifurcation  along  with  an 
interesting  sequence  of  period  doubling  bifurcation 
which  leads  to  chaotic  behavior.  These  issues  also  need 
further  investigation. 
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Fig.  1:  Stability  boundary  for  the  trivial  solution  (for 
P  =  -  0.5,  X  =  0.1  and  gb  =  0.5).  The  shaded  part 
indicates  the  stable  region,  ft,  1)  and  C2  indicate 
the  simple,  Hopf  and  Bogdanov  Takens 
bifurcation  varieties  respectively. 


Fig.  2:  Existence  of  the  non-trivial  steady  state  solutions 
for  p  =  -  0.5,  X  m  0.1,  gb  =  0.5  and  (a)  C2  *  1.,  d2  *  1. 
(b)  02  =  -  1.,  dj  =  1.,  The  number  inside  the  small 
ellipse  indicates  the  possible  number  of  non¬ 
trivial  equilibrium  solutions  that  may  exist  in 
that  region. 


F>g.  3:  Graph  of  Rp)  vs.  p.  Various  possibilities  are 
shown  for  the  existence  of  non-trivial 
equilibrium  solutions  of  r  for  the  nondegenerate 
case  (i.e.,  a  >  0).  (1).  b  <  0,  c  >  0,  D  >  0,  (2).  b  <  0, 
c<0,  (3).b>0,c<0,  (4). b<0, c>0, D<0. 

(The  dot  indicates  the  possible  non-trivial  steady 
state  solution.) 
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Fig.  4:  Curve  showing  various  branches  of  the  non¬ 
trivial  equilibrium  solution  for  different 
combinations  of  b  and  c,  as  i>2  is  varied  from  -1.0 
to  1.0. 


B  ig.  5:  Global  bifurcation  set  corresponding  to  Bogdanov 
Takens  bifurcation  variety. 
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Fig.  6:  Stability  boundary  for  the  trivial  solution  for 
P  *  -0.005,  X  =  0.0  and  gb  =  0.5.  The  shaded  pai  l 
indicates  the  stable  region.  The  Hopf  variety  <fi>) 
leads  to  the  Double  Hopf  bifurcation  variety  (J?2> 

as  p-+0. 
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*"'6  8:  Phase  plane  trajectories  (xj-  yj)  for  the  selected 
orbits  of  Fig.  6.  The  labels  in  this  diagram 
correspond  to  their  respective  positions  in  the 
bifurcation  diagram  (Fig.  6). 


Fig.  9:  Phase  plane  trajectories  [xy  yi)  for  various  values  of  u?;  (3  =  -0.5,  X  =  0.1, 
gb  *  0.5,  vi  =  0.1,  (a)  V2  -  0.3589,  (b)u2  =  0.35882,  (c)  v2  =  0.3588, 
(d)  V2  =  0.358785,  (e)  v?  =  0.35875,  (f)  Same  as  (e)  but  in  (x2-yi-y2)  plane, 
(g)  U2  =  0.3586,  (h)  vz  =  0.35855 
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Fig.  10:  (a).  Response  yj  as  a  function  of  u2;  P  =  *  0.1, 

A.  *  0.1,  gb  =  0.5,  \>i  =  -  0.3.  (b)  Same  as  (a)  but 
near  the  parameter  region  where  the  period 
doubling  occurs. 
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Fig.  11:  Phase  plane  trajectories  (xi  -  yj)  for  various 
values  of  t>2  for  the  upper  periodic  branch  of 
Fig.  10.  (P  =  -0.1,  A.  =  0.1,  gb  *  0.5,  ui  *  -  0.3), 
(a)  V2  =  0.357,  (b)  v2  «  0.355,  (c)  t>2  «  0.353, 
(d)  U2  =  0.352,  (e)  Same  as  (d)  but  in  iy\  -  y%-  *2) 
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The  effect  of  external  random  excitation  on  nonlinear  continuous  time  systems  is 
examined  using  the  concept  of  the  Lyapunov  exponent.  The  Lyapunov  exponent 
may  be  regarded  as  the  nonlinear/stochastic  analog  of  the  poles  of  a  linear  deter~ 
ministic  system.  It  is  shown  that  while  the  stationary  probability  density  function 
of  the  response  undergoes  qualitative  changes  ( bifurcations )  as  system  parameters 
are  varied,  these  bifurcations  are  not  reflected  by  changes  in  the  sign  of  the  Lyapunov 
exponent.  This  finding  does  not  support  recent  proposals  that  the  Lyapunov  exponent 
be  used  as  a  basis  for  a  rigorous  theory  of  stochastic  bifurcation. 


Introduction 

The  design  of  practical  control  systems  requires  that  the 
mathematical  model  used  be  sufficiently  robust.  This  implies 
that  qualitative  properties  should  be  preserved  under  the  effect 
of  all  possible  perturbations.  Such  systems  were  referred  to  by 
Andronov,  Vitt,  and  Chaikin  (1965)  as  “coarse  systems”  and 
formed  the  basis  for  the  concept  of  structural  stability  in  the 
mathematical  theory  of  dynamical  systems  (Arnold,  1983).  in 
the  context  of  single-input  single  o-itput  (SISO)  linear  control 
systems,  robustness  is  characterized  by  the  Nyquist  plot  and 
the  concepts  of  phase  and  gain  margins.  The  extension  of  these 
ideas  to  multiple-input  multiple-out  (MIMO)  linear  systems  is 
documented  by  Dorato  (1987).  Currently,  there  are  two  major 
approaches  to  robustness  analysis.  The  first  assumes  unstruc¬ 
tured  additive  or  multiplicative  perturbations  of  the  plant 
transfer  function.  Robustness  is  then  measured  by  the  singular 
values  of  the  return  difference  matrix.  This  is  a  frequency 
domain  technique  which  does  not  make  explicit  use  of  the 
nonlinearities  present  in  the  system.  The  second  method  is 
referred  to  as  a  structured  or  parametric  approach.  The  precise 
values  of  system  parameters  are  unknown  but  the  uncertainty 
is  assumed  to  be  bounded.  The  degree  of  robustness,  indicated 
by  the  system  poles  of  the  linearized  system,  is  then  determined 
in  terms  of  these  bounds  using  Kharitonov-type  theorems  (Bar- 
mish,  1988).  These  two  techniques  constitute  a  deterministic 
approach  to  the  robustness  analysis  of  linear  systems.  The 
possibility  of  random  parametric  or  external  excitation  is  not 
considered. 
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So  far  the  robustness  analysis  of  stochastic  systems  has  re¬ 
ceived  relatively  limited  attention.  Wonham  (1967)  investigated 
the  problem  of  a  linear  continuous  time  system  with  linear 
multiplicative  white  noise  gain  in  terms  of  optimal  stationary 
control  using  a  state  space  formulation.  A  similar  problem 
was  studied  by  Willems  and  Blankenship  (1971)  in  the  fre¬ 
quency  domain  (sec  also  Willems  and  Willems,  1983).  t.i  both 
instances,  robustness  was  measured  by  a  mean-square  type 
criterion.  Robustness  measures  of  stable,  linear  discrete  time 
systems  were  obtained  by  Yaz  and  Yildizbayrak  (1985)  and 
Yaz  (1988).  Their  robustness  measures  are  based  on  the  sample 
stability  of  the  perturbed  system.  The  references  mentioned 
above  are  mainly  concerned  with  the  effect  of  parametric  fluc¬ 
tuations/uncertainties  on  linear  systems.  For  nonlinear  sys¬ 
tems,  as  systems  parameters  are  varied,  qualitative  changes 
(bifurcations)  in  the  system  response  can  occur.  In  this  paper, 
the  robustness  of  such  bifurcations  under  the  effect  of  random 
external  excitation  will  be  examined.  In  the  following  section, 
the  problem  is  formulated  in  a  general  context  and  the  idea 
of  the  Lyapunov  exponent  as  a  quantitative  measure  of  ro¬ 
bustness  is  introduced.  Examples  of  systems  undergoing  co¬ 
dimension  1  and  2  bifurcations  are  then  studied. 


Problem  Formulation 

Consider  the  following  system: 

™  =  AX+f(X)  +  <n,U)  (1) 

where  A"  is  a  vector  in  R"  and  without  loss  of  generality  it  is 
assumed  that  A  is  in  canonical  form, /(AT)  are  nonlinear  terms, 
and  ij(f)  represents  independent  zero-mean  white  noise  exci¬ 
tation  of  unit  intensity.  Since  only  additive  noise  is  present,  it 
does  not  matter  whether  the  ltd  or  Stratonovich  interpretation 
is  used.  Assuming  that  a  stationary  solution  Xs  exists,  the 
problem  of  robustness  is  concerned  with  the  sample  stability 
of  this  stationary  solution.  As  with  deterministic  systems,  con- 
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sider  a  small  perturbation  x,  about  this  stationary  solution. 
Then  to  first  order  in  the  perturbation,  a  linear  system  with 
stochastic  coefficients  is  obtained: 


dx 

dt 


A  + 


V 

ax 


x=[A+F(t)]x. 


(2) 


Define  a  norm  p  =  (xTjr}'  \  the  exponential  growth  rate  of 
the  perturbation  is  then  given  by 

rf(logp)  xr[(/l  +  )  +  (A  +  F(t))T\x 

dt  ~  2p2  ‘  <3> 


Fig.  1  Lyapunov  exponent  (aott  tots  of  stabilify) 


The  Lyapunov  exponent,  A,  is  simply  the  time-averaged  ex¬ 
ponential  growth  rate  defined  by 


where 


/mo)  -  ('  jAM+fun  +  M+FW)1]* 

\p(0)/  J0  2p: 


(4) 


ds. 


The  stationary  solution  is  sample  stable  if  A  <  0,  and  this  is 
a  necessary  condition  for  a  robust  system.  Assuming  that  the 
stationary  state  is  also  ergodic,  the  temporal  average  may  be 
replaced  by  the  ensemble  average  and  a  quantitative  measure 
for  robustness  can  then  be  computed,  i.e., 

VU;4+F(f))+M+f(0)V 


■** 


2p2 


x 


(5) 


One  may  regard  the  Lyapunov  exponent  as  the  stochastic  an¬ 
alog  of  the  poles  of  a  deterministic  system.  This  concept  of 
the  Lyapunov  exponent  (Bylov  et  al.,  1966)  has  been  proposed 
by  L.  Arnold  (1988)  as  a  basis  for  a  rigorous  theory  of  sto¬ 
chastic  bifurcation.  A  similar  approach  was  also  used  by 
Caughey  and  Gray  (1965),  Infante  ( 1968),  Kozin  and  Wu  (1973), 
and  Ariaratnam  and  Xie  (1988a)  in  their  derivation  of  sufficient 
conditions  for  the  sample  stability  of  linear  systems  with  sto¬ 
chastic  coefficients. 

One-Dimensional  Systems 

For  one-dimensional  systems,  the  computation  of  the  Lya¬ 
punov  exponent  is  fairly  simple.  Consider  the  general  one¬ 
dimensional  system  written  as 

~^^\X+f(X)  +<m(r),  (6) 


and  the  perturbed  system  is 


X  + 


df(Xs) 

dX 


x. 


(7) 


I  f  the  stationary  state  is  ergodic,  then  the  Lyapunov  exponent 
is  given  by 


A  =  X+£ 


df{Xs)  1 

.  dX  J 


(8) 


and  the  stationary  probability  density  function  can  be  obtained 
by  solving  the  Fokker-Planck  equation  (FPE)  for  the  system 

9  \l)X+AX)]p(X)-^^[.  (9) 


0= 


ax 


b 


The  stationary  probability  density  function  (pdf)  is  given  by 

pt(X)=Ntxp{p  ^Xy+j  fiz)dzj'j  (10) 

where  N  is  an  appropriate  normalization  constant.  It  is  in¬ 
structive  to  consider  a  simple  example: 

Example  1  (Simple  Bifurcation— Soft  Loss  of  Stability). 
For  a  symmetric  system  undergoing  a  soft  loss  of  stability  (i.e.. 


P*(*> 


Ftp.  2(4  Probability  density  function  (soft  lost  of  stability),  *  =  -  0.2 


Ps(*> 


Fig.  2 (4  Probability  density  function  (soft  loss  ol  stability),  *  =  02 


a  loss  of  stability  to  a  nearby  equilibrium,  f(X)  =  aXx  with 
a  <  0)  and  the  stationary  probability  density  function  is 


p,(X )  =  exp  ^  (AA^/2  +  aX*/ 4)^  j 


j  exp  (p  (XX2 /2  +  aX*/ 4)^  dX 


(ID 


Assuming  that  the  stationary  solution  is  ergodic,  the  Lya 
punov  exponent  is  then  given  by  Eq.  (8): 


A  =  A  +  3a£-[A'‘l. 


(12) 


The  Lyapunov  exponent  A  is  plotted  against  the  eigenvalue 
X  for  the  deterministic  system  in  Fig.  1  with  <r  =  0. 1  and  a 
=  -  I.  The  deterministic  system  undergoes  a  bifurcation  at  X 
=  0  and  the  probability  density  function  changes  from  a  un- 
imodal  density  to  a  bimodal  density  (Fig.  2(o)  and  2(6))  As 
evident  from  Fig.  I ,  the  Lyapunov  exponent  A  is  negative  u  hich 
indicates  that  this  is  a  robust  feature,  but  the  system  is  least 
robust  at  X  =  0.  Physically,  this  implies  that  perturbations 
will  take  a  much  longer  time  to  decay.  It  should  be  noted  ihat 
the  Lyapunov  exponent  does  not  indicate  the  qualitative  change 
in  the  stationary  pdf.  This  point  was  not  adequately  empha¬ 
sized  by  L.  Arnold  (1988)  in  his  theory  of  stochastic  bifurcation 
based  on  the  Lyapunov  exponent. 
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Fig.  3  Lyapunov  exponent  (hard  tote  of  stability) 

Ps!*) 


Fig.  4(e)  Probabilily  density  (unction  (hard  lots  ol  stability),  X  =  - 1.5 


ps;*) 


Fig.  4(b)  Probability  density  (unction  (hard  loss  ol  Stability),  X  =  -  0.9 


Example  2  (Simple  Bifurcation — Hard  Loss  of  Stability). 

Now  consider  the  one-dimensionai  system 


dX 

dt 


=  xjt'  +  cyr-  +  cs*s+<ji)<o 


(13) 


«  here  C,  >  0,  C<  <  0  for  a  hard  loss  of  stability.  The  steady- 
stale  pdf  is  given  by 


2  /  X2  X4  „  X*\~ 

exp[P  (xT+C3T+C5t) 


PAX )  =  A/  exp  j 
where  the  normalization  constant  N  is  defined  by 
N= 


(14) 


f®  2 

.  (  X2  „  X*  „  x\ 

1  exp  - 

J .  „  0 

2  (xT+C}T+CjT)j 

dX  (15) 


and  the  Lyapunov  exponent  is  (from  Eq.  (8)) 

A  =  X  +  3Cj£1*3)  +  5C$£lJf*]  ( 1 6) 

where  the  expectation  is  taken  with  respect  to  the  steady-state 
pdf  defined  by  Eq.  (14).  Numerical  results  for  C3  *  2,  C5  = 
-  1  and  a1  ~  1  are  shown  as  Fig.  3.  In  this  case,  the  system 
is  least  robust  just  before  the  deterministic  bifurcation  point 
at  X  =  0.  It  is  interesting  to  compare  this  result  with  the 


Fig.  S  Lyapunov  axponant  tor  dX  =  (1  ♦  XX  -  XVH  ♦  s/iXo  dW 


qualitative  changes  in  the  pdf  as  X  is  varied  (Fig.  4).  It  can  be 
seen  that  the  decrease  in  robustness  coincides  with  the  emerg¬ 
ence  of  two  additional  peaks  in  the  pdf.  Once  again,  it  should 
be  noted  that  the  Lyapunov  exponent  remains  negative  and 
cannot  be  used  as  a  bifurcation  parameter. 

The  results  derived  here  for  nonlinear  systems  under  external 
random  excitation  is  by  no  means  typical.  In  fact,  parametically 
excited  systems  studied  by  previous  researchers  in  stochastic 
bifurcation  theory  are  extremely  sensitive  to  errors  in  the  ref¬ 
erence  input.  Consider  the  nonlinear  system  with  fluctuating 
gain: 

dX=(\X-X})dt  +  \/2XodW,  (17) 

where  "o”  denotes  that  the  Stratonovich  interpretation  is 
used.  It  is  well  known  that  the  Lyapunov  exponent  is  simply 
X  and  the  trivial  solution  X  =  0  loses  its  stability  at  X  =  0 
and  this  coincides  with  a  qualitative  change  in  the  probability 
density  function.  Suppose  now  that  a  nonzero  reference  input 
is  present,  i.e., 

dX^V  +  XX-X^dt  +  y/lXodW.  (18) 

Since  the  drift  term  is  not  zero  at  X  =  0,  the  system  is  ergodic 
and  a  steady-state  solution  exists  and  is  defined  by  the  steady- 
state  pdf: 

(X*-  ‘exp[  ~(\/X  +  X2/2)]/C  X>0 
PAX)  =  (19) 

LO  else 

where 

C=  (  A* "  'exp[  -(1/X  +  Xl/2))dx, 

Perturbing  about  this  stationary  (ergodic)  solution  leads  to 
the  following  linear  system  with  stochastic  coefficients: 

<tr  =  (A  -  l(X,)2)xd!  +  sfixodW.  (20) 

The  Lyapunov  exponent  is  then  given  by 

A^?!l0«7$rr-x-3£,jr’1  <2I> 

since  the  Wiener  process  W(l)  -  (/loglogO'^asf  —  »  with 
probability  one.  The  Lyapunov  exponent  Eq.  (21)  is  plotted 
in  Fig.  5  and  is  negative.  The  variation  of  the  Lyapunov  ex¬ 
ponent  with  X  is  not  unlike  that  of  a  system  perturbed  by 
external  random  excitation.  This  “quasi-external  excitation” 
effect  of  a  nonzero  reference  input  is  apparent  if  the  trans¬ 
formation  X  =  Y  +  C  is  made  in  Eq.  (18)  so  that  the  nonzero 
reference  input  is  eliminated.  One  sees  that  an  additional  ex¬ 
ternal  excitation  term  (CdW)  is  then  generated  from  the  mul¬ 
tiplicative  noise  term  XodW. 


Two-Dimensional  Systems 

Noting  that  the  problem  of  robustness  leads  to  a  linear 
system  with  stochastic  coefficients,  it  is  convenient  to  first 
consider  the  general  two-dimensional  linear  system  with  sto¬ 
chastic  coefficients: 
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“  =  Fh  ( # )X|  +Fi2{t)X2 


■=\n  +  2a(X1+Yi) 


-jjf  =  F2 1  ( t)x ,  +  Fn  ( t)x2 

where  the  Fu(t)  are  colored  noise  sources  generated  from  white 
noise.  Making  the  transformation,  xt  «  rcos0,  x2  *  r  sine, 

j  =  \  l(F„  +  ^22)  +  (Fl2  +  Fu)sin23  +  (F„ -Fu)cos2B)r 

-^*2  K^2i  “^12)  +  <^22  —  )s»o2tf  +  (Fn  +F|2)cos2tfJ  (23) 

If  the  0  process  and  the  Fv(t)’s  are  ergodic  then  the  Lyapunov 
exponent  is  given  by 

A  =  -  E[(Fu  +  fa)  +  yj  (fi2  +  F2  i  )J  +  (f ii  -  fa)Jsin(20  +  0)] 


where  4  is  a  random  phase  angle  defined  by 


(F»-Fa\ 
\FI2  +  F2l) 


<i>  =  arctan 


Noting  that  lsin(20+$)l  ss  1,  an  upper  bound  for  the  Lya¬ 
punov  exponent  can  be  found  and  hence  a  sufficient  condition 
for  sample  stability  is: 

£l(fn  +fa)  +-\/  (F^2  +  F2l)2  +  (fu-fal^cO.  (26) 
A  necessary  condition  for  sample  stability  requires  that  the 
joint  moment  of  the  stochastic  coefficients  and  the  phase  proc¬ 
ess  0  be  computed.  These  ideas  are  applicable  to  the  following 
system: 

Example  3:  Hopf  Bifurcation  (Dynamic  Instability). 
Consider  the  nonlinear  system  perturbed  by  independent  ex¬ 
ternal  white  noise  sources: 

^=nX-u>Y+  (aX-bY)(X2  +  yJ)  +  <nj|(0 


^-  =  u>X+,iY+(aY+bX)(X2+Y2)+<ty,2(t)  (27) 

at 

where  the  deterministic  terms  correspond  to  the  normal  form 
for  the  Hopf  bifurcation.  A  stationary  state  exists  for  the  case 
a  <  0,  and  the  pdf  is  given  by 

<»> 

and 

,  r  f"  T2  (  (Jf'+y2) 

2 

+  a  ■*  jY  )ZSj  j dXdY.  (29) 

Perturbing  about  this  stationary  state,  the  resulting  linear  sys¬ 
tem,  Eq.  (22),  is  defined  by  x,  =  x,  x2  =  y  and 

Flt  =  n  +  aOX2+  Y2)-2bXY 

-u>-b(X1+iY1)  +  laXY 

F2l  =  u>  +  b(3X2  +  Y1)  +  2aXY 

Fn  =  ii  +  a(X*  +  2Y*)  +2bXY  (30) 

which  leads  to 


+  yJ<t  +  b2{(X2-  yJ)sintf-2ATeosC!!r  (Ji, 

—  =  v  +  2b(X1+Yl)  +  \Ja2  +  b2[(X2-  Y:) cov;  -  2ATsin-;i 

where  ^  =  2 1  +  arcta n(a/b).  Since  w  is  nonzero  and  the  X 
Y processes  are  assumed  to  be  ergodic,  the  Lyapunov  exponent 
is  then  given  by 

A  =  Eht  +  2a(X2+  Y1) 


+  \/a1  +  b2\(X2- V2)sin^  -  ZATVcosCl] 


where  the  expectation  is  taken  with  respect  to  the  steady-state 
joint  probability  density  function  ps( X,  T.  6).  A  perturbation 
expansion  for  this  steady-state  pdf  may  be  conducted  usme 
a  separation  of  time  scales  technique  given  in  Blankenship  and 
Papanicoloau  (1978).  Roughly,  since  the  X,  Y  processes  are 
assumed  to  be  stationary  to  start  with,  these  processes  must 
have  evolved  sufficiently  fast  to  reach  a  stationary  state  and 
should  be  scaled  accordingly  for  consistency.  The  technique 
has  also  been  applied  by  Horsihemke  and  Lefever  (1984)  and 
is  explained  in  Appendix  A.  Then  the  Lyapunov  exponent  can 
be  computed  term-wise  from  such  an  expansion  has  been  pros  en 
by  Arnold  et  al.  (1988).  For  simplicity,  the  technique  will  be 
applied  for  the  case  b  -  0  and  the  angular  frequency  u;  is 
large.  The  system  defining  the  steady-state  pdf  is  scaled  as 
follows: 

l(*:-  T2)sin20  -  2XYco$20] 

at  €  € 


biX-wY+aX(X2+Y2)]  +  ~  *,(/) 


~=4l wX-»Y+aY(X2+Y2))+-V2U)  (33) 

at  t  t 

Let  the  steady-state  pdf  take  the  following  form: 

PAX,  Y,  8)=po(X,  Y,  6)  +  tpAX.  Y,  8)+ -  (34) 

The  steady-state  Fokker-Planck  Equation  for  Eq.  (33)  may  be 
written  as 

p4olp,l  +  ^Z.  .1*1=0,  <3?> 

where  the  operators  Lo  and  L,  are  defined  as  follows: 

ioW-  ([M*-wy+a*(*J+  y2)l lf>-y  f|) 

~dY  [[^  +  ^Y+aY(X2+Y‘)]p-1^j  -  * 

Iilp]=  -a4z  [\(X2-  yJ)sin20 -  2ATcos201p).  (-’<•> 

ou 

Substituting  Eq.  (34)  for  the  pdf  leads  to  the  system  ot 
equations: 

LolPol  =  0 

LoIp,l  =  LiU»o] 

«  •  «  • 

LolP,l  =  ,J” 

The  solution  to  the  lowest  order  equation  is 
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PAX,  V.  6)=PAX,  Y)pAB) 

(A ,2+Yi) 
l - - - 4  a 


«  N  exp 


{?(' 


a',  ry\l  j_ 

4  )\lt' 


(38) 


and  an  approximation  to  the  Lyapunov  exponent  is  then  ob¬ 
tained  from  Eq.  (32)  by  taking  expectation  with  respect  to  Eq. 
(38), 

A=*  +  2a£lA':+  Y*].  (39) 


This  is  compared  in  Fig.  6  with  the  sufficient  condition  ob¬ 
tained  from  Eq.  (26),  which  states  that  the  system  is  robust  if 

H  +  (2a  +  la  I  )£IX2  +  Y1]<0.  (40) 

For  a  =  -0.5  and  a2  =  I,  the  sufficient  condition  Eq.  (40) 
is  satisfied.  This  indicates  that  the  system  is  robust  but  as 
expected,  it  provides  a  more  conservative  robustness  measure 
compared  to  Eq.  (39). 

One  may  also  check  if  the  polar  representation  of  the  Hopf 
bifurcation,  Eq.  (27),  is  robust  to  external  excitation.  Consider 


Fig.  •  Lyapunov  exponent  (Hop!  bifurcation)  (1)  uppor  bound— non- 
radlal  axcftatlon,  (3)  axact— radial  aacltatfon,  (X)  lowaal-ordar  approxi¬ 
mation—  nonradial  axcltatton 


dR 

dt 


=  m*  +  c/?J  +  eRZ 2  +  x/2,,  (/) 


—  =  n/(  +  d/?3  +  diji(r),  a<0, /faO 

dt 


~  =  \Z+bZ)  +  dR2Z  +  y/2  vAO, 


(45) 


—  =  u  + bR2  +  or)2(t)  (41) 

at 

where  the  noise  terms  are  included  for  the  robustness  analysis 
and,  for  purposes  of  comparison,  b  will  be  set  to  zero.  Ob¬ 
serving  that  the  amplitude  R  is  decoupled  from  the  phase,  the 
Lyapunov  exponent  is  then  given  by 

A  =  M  +  3a£Tf?J],  (42) 


where  the  expectation  is  taken  with  respect  to  the  probability 
density  function 

PAR.  9)  =  Nexp0  (/tR2/2+aR*/4)j  — 


and 


AT 


1  =  |  exp  0  (»R:/2  +  aR*/4)j dR. 


(43) 


The  Lyapunov  exponent,  Eq.  (42),  is  also  plotted  against  p  in 
Fig.  6  for  the  same  parameters,  i.e.,  a2  =  1  and  a  =  -0.5. 
It  is  evident  that  the  polar  representation  is  less  robust  to 
external  random  excitation.  Physically,  in  Eq.  (41),  the  per¬ 
turbations  are  aligned  with  the  radial  direction  and  hence  the 
perturbations  should  have  a  greater  effect  on  the  amplitude. 

Example  4  (Hopf-Pitchfork  Interaction— Coupled  Dy- 
namic/Stafic  Instability).  Consider  the  deterministic  normal 
form  for  the  Hopf-pitchfork  bifurcation  perturbed  by  external 
independent  white  noise  sources: 

=  nR  +  c '  R}  +  e '  RZ2  +  «|i),  ( t ) 


^^-kZ+b'Z'  +  dRtZ  +  OtfAl),  (44) 

dt 

where  R  represents  the  dynamic  mode  (R  >  0),  and  Z,  the 
static  mode.  Based  on  Example  3,  the  polar  representation  is 
used  since  from  a  robustness  viewpoint,  this  corresponds  to  a 
"worst-case”  situation.  It  is  convenient  to  rescale  the  state 
variables  by 


R—j2R/o,,  Z-yfiZ/ot, 

resulting  in 


where  the  new  system  parameters  are  given  by 
h=6'<<rj)J/2,  c  =  c'(o,)1/2,rf*<f'(o,)1/2.e«e'(ff2),/2. 

An  explicit  normalizable  solution  to  the  steady-state  FPE 
for  Eq.  (45)  can  be  found  for  the  case  d  =  e  =  k,  provided 
b.c  <  0  and  if  k  >  0,  then  be  -  k2  >  0  must  also  be  satisfied. 
The  steady-state  probability  density  function  is  then 

PAR,  Z)  =  Ncx p0  (2p  +  cR2)R2  +  2kR*Z2  +  (2X  +  bZ2)Z2j, 

(46) 

with  S  as  a  normalization  constant.  Perturbing  about  this 
stationary  state  and  to  first  order  in  the  perturbation, 

(M  +  3cfl2  +  *Z2)r+  (2*J?Z)z 


(2**Z)r+  (\  +  IbZ2  +  kR2)z.  (47) 

at 

Using  the  formula  derived  earlier,  Eqs.  (22M23),  and  letting 
r  =  pcostf,  z  =  psintf,  the  growth  rate  of  the  norm  of  the 
perturbation,  p,  is  governed  by 

[(X  +  p)  +  (3c +  *)**+  (3f>  +  *)Z,  +  4*j?Zsin2d 

at  2 

+  (0»-X)  +  (3c-*)/fI+(*-3h)Z,)cos2dlp.  (48) 
Scaling  the  9,  R,  Z  processes  as  in  Appendix  A  yields 

?  =  ^  (X  -  p)sin2 9  +  r  1 4kRZcos2$ 
at  2  2e 

+  [(3  b  -k)Z1+(k~  3c)R2)$in29  J 

^  -T  (m*  +  c*J  +  **Z2 )  +  -  >/ 2  „(f) 
at  t  t 

*=  T  (XZ  +  bZ1  +  kR2Z)  +  -  V2^0-  (49) 

at  t  t 

For  simplicity,  let  *  ■  3ft  *  3c.  The  steady-state  FPE  is 
given  by 

(50) 

with  Lo,  Lu  Li  defined  by 


4 


i 


I 

1 
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£,lpl  =  2*/?Z~  [pcos28] 
qQ 

L:[pl*“"l/»in28). 

Let  the  perturbation  expansion  be  of  the  form 
p,(R,  Z,  8)  =*,(*.  2,  8)  +  ep,  (/?,  Z,  B) 

+  <WR.2,e)+...,  (31) 

and  the  system  of  perturbation  equation  is 
LolPt>]  =  0 

Lo(Pil  =  L|(pol 
Lolpd  -  Z-  ilpil  +  /-2lA»l 


P(R,  Z)  -  -  M(Q  2)  +  M(2  0) 


drm 

3 

r  3f,i 

P 5  dR 

3Z 

[p’dz\ 

which  is  a  linear  self-adjoint  operator.  Substituting  into  the 
FPE,  Eq.  (5)  and  collecting  terms,  the  0(<°)  equation  is  trivally 
satisfied  with  a,(«.Z.»)  =  p,(«,  Z)r0(  8),  but  unlike  Example 
3,  r<£9)  has  yet  to  be  defined.  The  (Xt1)  equation  is 

LolP/il  =  IiUYo(8)]  =  2 kRZp,(R,  Z)  ^  =  ro(8)cos28.  (54) 

00 

According  to  the  solvability  condition, 
j  {  J  2JcRZp,(R,  Z)  ^  =  r„(8)cos28 

xpAR,  Z)N(B)dZdRd9  =  0 

where  the  term  p,(R,  Z)S(9)  belongs  to  the  kernel  of  L0. 
Since  this  has  to  hold  for  an  arbitrary  N(6),  it  is  required  that 


rr. 


2JcRZp]{R,  Z)dZdR  =  0. 


As  p,(R,  Z)  is  symmetric  in  Z,  this  implies  that  Fredholm’s 
alternative  is  satisfied  and  that  a  unique  solution  for  r,  can  be 
found.  Using  the  fact  that  L0  is  a  linear  operator,  r,  can  then 
be  expressed  as  the  sum  of  a  homogeneous  and  particular 
solution 

r,(K,  Z,  e)  =  H($)  +  2k  ^=r„(8)cos2M>(*,  Z),  (56) 

dv 

where  H(9)  is  the  homogeneous  solution  of  L0[pjf(6) I  =  0 
and  P(R,  Z)  is  a  particular  solution  of  Lo\pf{R,  Z)]  = 
RZp,(R,  Z).  It  remains  now  to  solve  for  P(R,  Z).  Using 
Galerkin’s  approximation,  let  P  =  £CmmRmZ'  and  take  weight 
functions  of  the  form  RPZ‘.  Define 

A/<«,«)=(  f  RmrPs(R,  Z)dZdR.  (57) 

A  simple  3-mode  expansion,  P(R,  Z)  =  CjpR2  +  CURZ 
+  GnZ\  yields 


Now,  at  0(el): 


(X  -  n)  d 


Lol/Vi)  =  2 kRZ  —  [/-,(/?.  Z,  8)cos28|  +  — —  —  (/-SJ<e)sin2&) . 


Substituting  for  r,(/?.  Z,  B)  and  invoking  Fredholm's  alter¬ 
native,  r^B)  has  to  satisfy 


Lob»l  =  L,lp«-i)+LilA.-2].  (52) 

It  is  convenient  to  let  p„(R,  Z,  0 )  =  PAR.  Z)r„(r,  Z,  0), 
where  ps(R,  Z)  is  defined  in  Eq.  (46),  so  that  the  operator 
La\p.\  takes  the  form 


—  pOo(8)sin28  +  cos20  —  (r„<ekos20]  j  =0,  (60) 

with  C  =  (X  -  n)/(\6lcE\RZP{R,  Z».  where  the  expectation 
is  taken  with  respect  to  ps(R,  Z).  This  equation  may  be  re¬ 
written  as 

3  f  r  i  i  )  i  d~- 
~dB  ~  Cs,n2fl  +  2 S,n<W  +2  3?'V°S'2(l'i0'  (61 1 

which  takes  the  form  of  a  steady-state  Fokker-Planck  equation 
with  drift  and  diffusion  defined,  respectively,  by 

+(8)  =  -  £  Csin20  +  ^  sin48  and  ♦:(8)  =  cos:28.  (62) 

The  process  evolves  on  the  half  circle  8  €  [0,  x]  with  sin¬ 
gularities  at  8  =  x/4  and  3x/4.  As  C  *  0,  the  drift  at  the 
singular  points  is  not  zero  and  thus  the  singular  points  can  be 
classified  as  left  and  right  shunts  depending  on  the  sign  of  the 
drift  term.  The  boundary  point  classification  is  performed 
using  the  Feller  classification  scheme  (1954).  It  is  shown  in 
Appendix  B  that  the  boundary  points  8  =  0,  x  are  both  per¬ 
fectly  reflecting  regular  points.  Therefore,  the  8  process  is 
ergodic  and  the  Lyapunov  exponent  can  be  found  using  Eq 
(24)  to  be 

A  =  — ~  +  £{cos28]  +  *(£■{/?•]  +  £[Z:j).  (63) 

Now  consider  the  case  C  >  0,  (similar  results  can  be  derived 
for  C  <  0).  The  singular  points  are  given  as  left  and  right 
shunts  for  8  =  x/4  and  3x/4,  respectively.  That  a  nontrivial 
normalizable  solution  for  r0  exists  on  the  intervals,  (0,  x/4) 
and  (3t/4,  *)  is  shown  in  Appendix  B.  r„  is  given  by: 

C  >  0. 

r  n 

- — - — — —  0  <  8  S  *74  or  3x/4  £  8  S  * 

r0(8)  =  exp(C/cos28)cos28 

0  else 

V  (64) 

where  Afis  a  normalization  constant.  The  Lyapunov  exponent, 
Eq.  (63),  is  plotted  against  *  (where  X  =  v  and  fi  -  v  +  1) 
in  Fig.  7  for  k  -  -I.  The  deterministic  system  undergoes 
primary  bifurcations  at  ¥  -  Oand  -  1.  A  secondary  bifurcation 
occurs  at  r  =  0.5.  It  can  be  seen  that  these  features  are  robust 
(i.e.,  the  Lyapunov  exponent  remains  negative)  and  a  decrease 
in  robustness  coincides  with  the  occurrence  of  the  first  prir  ry 


Fig.  7  Lyapunov  axponant  (Hop!  pticMof*  bifurcation) 
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bifurcation.  Once  again.  the  Lyapunov  exponent  does  not  re¬ 
flect  the  qualitative  changes  in  the  pdf  but  does  indicate  when 
nonlinear  effects  become  significant. 


Conclusion 

The  major  conclusions  of  this  investigation  are  as  follows: 

1  While  the  Lyapunov  exponent  provides  a  robustness  meas¬ 
ure  for  the  combined  effects  of  nonlinearities  and  noise,  it  is 
not  a  suitable  bifurcation  parameter  for  nonlinear  stochastic 
systems  in  the  sense  that  it  does  not  reflect  the  qualitative 
changes  in  the  steady-state  probability  density  function.  This 
occurs  whenever  X  =  0  is  nor  a  solution  of  the  nonlinear 
stochastic  system.  Examples  of  nonlinear  systems  under  ex¬ 
ternal  (additive)  random  excitation  and  systems  under  para¬ 
metric  random  excitation  (multiplicative  noise)  with  a  nonzero 
reference  input  were  shown  to  possess  negative  Lyapunov  ex¬ 
ponents. 

2  The  results  of  this  investigation  suggest  that  the  extremum 
of  the  Lyapunov  exponent  (with  respect  to  the  bifurcation 
parameter  for  the  deterministic  system)  provides  an  indication 
of  when  nonlinear  effects  become  important.  For  the  one  and 
two-dimensional  systems  perturbed  by  external  random  exci¬ 
tation,  it  was  observed  that  this  extremum  follows  the  occur¬ 
rence  of  the  first  bifurcation  in  the  corresponding  deterministic 
system  It  is  left  as  a  conjecture  that  this  is  a  general  char¬ 
acteristic  of  nonlinear  systems  perturbed  by  external  random 
excitation. 
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APPENDIX  A 

The  scaling  employed  is  adapted  from  Horsthemke  and  Le¬ 
fever  (1984)  and  may  be  justified  as  follows: 

1  The  X  and  Y  processes  are  assumed  from  the  start  to  be 
stationary  (ergodic)  processes.  This  stationary  state  is  achieved 
only  in  the  limit  /  —  ot>.  Hence,  in  the  Ito  representation  of 
the  X  and  Y  processes,  time  is  scaled  as  I  —  l/t2  so  that  the 
stationary  state  is  reached  in  the  limit  e  —  0. 

2  Now  consider  the  effect  of  such  a  scaled  process  Z(t)  = 
X(i/e2)  in  the  8  equation.  For  example, 

d6=f(8)dt+Z(l)g(8)dt.  (65) 

The  spectrum  of  Z(t)  =  Xd/t1)  is  given  by 
l  r” 

Sz(<*i)  =  ~  j  E\Z{t)Z(t  +  T)]exp(-juT)dT 

\  £]*< t/t2)X{ l/t2  +  r/f2 )]exp(  - jtar)di 

=^j  Bvnt/t)xv/*+'n 

xexp( -yf2uis)<fc(let  s=r/e2) 
=  t2SxU2u)  (66) 

which  tends  to  zero  as  *  —  0.  This  implies  that  by  simply 
scaling  time  alone,  one  will  obtain  a  noiseless  limit.  Hence, 
the  amplitude  of  the  process  X(t/e2)  should  be  multiplied  by 
a  factor  of  1/e  to  ensure  that  noise  effects  are  not  excluded  in 
the  limit. 


APPENDIX  B 

The  ergodic  assumption  is  shown  to  be  valid  by  considering 
the  properties  of  the  sample  paths  of  the  8  process  on  the  unit 
half  circle.  The  drift  and  diffusion  are,  respectively, 

*(8)  =  -  Csin2#+  ^  sin-tfj  and  *2(8) = cosJ2*.  (67) 

Take  C  >  0,  and  by  symmetry  it  is  only  necessary  to  study 
the  process  on  (0,  jr/4).  The  diffusion  is  singular  at  I  *=  t/4 
with  #(r/4)  =  -C  <  0,  yielding  8  =  t/4  as  a  left  shunt. 

In  order  to  complete  the  proof  of  ergodirity,  the  behavior 
at  8  -  0  must  be  known.  To  this  end,  we  introduce  scale  and 
speed  measures  defined,  respectively,  as 

S(0)=j  s(x)dx,  A/(0)*j  m(x)dx  (68) 

where 


s(x)**Expl-J>(x)}, 


m{x) 


1 

P(x)s(x) 


(69) 
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s(jc)  and  m(x)  are  called  the  scale  and  speed  densities.  The 
scale  and  speed  measures  are  found  to  be 

E*P[cos2ij]  f,EXP[  cos2,»] 

»«-)  lcoU,l  "“"J  Ico, 2,1 

(70) 

Again,  only  consider  C  >  0  and  8  €  [0,  t/4).  Both  functions 
are  bounded  and  continuous  on  8  €  (0,  e/4)  and  as  such  the 
integrals  on  this  region  are  also  bounded.  Thus,  8  =  0  is  a 
regular  point. 

The  property  of  pure  reflection  can  be  determined  from  a 
physical  argument.  It  is  clear  that  unless  the  coupling  coeffi¬ 
cients  are  both  zero,  it  is  not  possible  to  have  a  solution  with 


?(/)  =  0.  Since  9  =  0  implies  z  =  0.  there  can  be  no  accu- 
mutation  of  probability  mass  at  the  point. 

As  8  =  0  is  purely  reflecting  and  8  =  t/4  is  a  left  shunt 
the  region  (0,  t/4)  has  the  property  of  recurrence  ot  ergodicity 
and  a  unique,  normalizable  stationary  density  exists  on  this 
interval.  Similar  conclusions  may  be  made  for  the  interval 
(3t/4,  t).  Should  8  start  outside  the  intervals  (0.  t/4)  and 
(3r/4,  t),  it  can  be  shown  (Nishoka,  1976)  that  the  8  process 
moves  into  either  one  of  the  two  intervals  in  finite  time  with 
probability  one.  The  temporal  average  can  then  be  replaced 
by  the  ensemble  average  computed  with  respect  to 

p(8)  =  Nm(6)  =  CQ%2eExp[c7coi2f)\  <7» 

where  N  is  a  normalization  constant.  Similar  results  can  be 
obtained  for  C  <  0. 


> 


1022  /  Vol.  59,  DECEMBER  1992 


Transactions  of  the  ASME 


